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CHAPTER  I 


INTRODUCTION 

Solving  large  scale  mathematical  problems  has  been  and  will  always  be  a 
vital  concern  in  many  areas  of  life.  This  concern  has  led  to  an  increased  popular¬ 
ity  and  explosive  growth  of  numerical  analysis.  Numerical  analysis  is  the  theory 
of  constructing  methods  for  approximating,  in  an  efficient  manner,  the  solutions 
to  mathematical  problems.  Existing  methods  are  grouped  into  two  types:  direct 
methods  and  iterative  methods  (Golub  and  Vein  Loan,  1983).  Direct  methods  will 
determine  a  solution  exactly,  up  to  the  precision  of  accuracy  of  the  computer,  in  a 
finite  number  of  steps.  The  goal  of  iterative  methods  is  to  start  with  an  initial  guess 
of  the  solution  and  then  improve  the  guess  by  the  use  of  an  updating  step  which 
is  called  an  iteration.  A  succession  of  iterations  will  produce  a  sequence  of  scalars 
-  real  or  complex  numbers  -  or  vectors,  as  appropriate,  which  converges  to  a  limit, 
the  solution.  If  the  limit  is  not  obtained  in  a  finite  number  of  iterations,  it  may  be 
approximated  to  a  desired  accuracy  after  a  finite  number  of  iterations.  This  study 
will  be  concerned  only  with  iterative  methods. 

Before  the  age  of  digital  computers,  methods  requiring  a  large  amount  of  com¬ 
putational  effort  were  impractical,  if  not  totally  unreasonable,  to  apply.  However, 
we  now  have  the  high  speed  computers  available  that  allow  us  to  solve  large  scale 
problems  to  a  high  degree  of  accuracy.  But  we  still  have  a  problem:  computer  time 
costs.  Therefore,  it  was  essential  that  the  subject  area  of  numerical  analysis  respond 
to  the  overwhelming  need  for  faster  computation  and  reduced  computer  time.  The 
result  was  the  development  of  sophisticated  numerical  methods  called  acceleration 
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techniques.  These  techniques  use  the  original  sequence  to  produce  a  new  sequence 
which  converges  to  the  same  limit  as  the  original  one,  only  faster. 

An  acceleration  technique  is  presented  in  the  form  of  what  is  called  an  ex¬ 
trapolation  algorithm  (Skelboe,  1980).  An  extrapolation  algorithm  determines  an 
element  of  the  second  sequence  from  some  desired  number,  say  k,  of  consecutive 
elements  of  the  original  sequence.  In  other  words,  assume  our  original  sequence  is 

{®n}  —  {xo,  ®l ,  X2,  .  .  .}. 

We  then  extrapolate  (compute)  a  second  sequence,  {y„},  by  applying  the  algorithm 
to  the  k  terms  and  represent  an  extrapolation  by  the  notation 

2/p+k-i  =  EXT(xp,  Xp+i,  Xp+2,  •  •  • ,  Xp+k- 1),  p  —  0,1,2,... 

The  purpose  of  this  study  is  to  compare  several  of  these  acceleration  tech¬ 
niques.  The  techniques  studied  will  include  Wynn’s  (1956)  epsilon  algorithm,  the 
modified  epsilon  algorithm  (Cheng  and  Hafez,  1959),  Cabay  and  Jackson’s  (1976) 
Minimal  Polynomial  Extrapolation  method  (MPE),  a  matrix  Full  Rank  Extrapola¬ 
tion  (FRE)  originally  developed  by  Henrici  (1964)  and  modified  to  a  Reduced  Rank 
Extrapolation  (RRE)  by  Eddy  (1979),  and  Anderson’s  (1965)  generalized  secant 
methods.  In  addition,  a  method  derived  by  Aitken  (1936-37)  for  scalars  and  mod¬ 
ified  for  vectors  by  Jennings  (1971)  will  be  studied  in  the  early  chapters  to  help 
establish  notation  and  to  set  the  pattern  of  how  extrapolation  algorithms  will  be 
developed. 

The  focus  of  this  study  will  be  on  acceleration  techniques  for  finding  the  so¬ 
lution  of  the  problem 

F(x)  =  0,  (1) 

where  F  is  an  operator  on  a  scalar  x  or  on  a  m-dimensional  vector  x.  The  theory 
will  be  developed  in  the  beginning  with  the  use  of  scalars;  however,  all  theory  will 
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quickly  be  related  to  vectors  and  by  the  end  of  the  study  our  only  concern  will  be 
solving  (1)  operating  on  m-dimensional  vectors. 

For  convenience,  Equation  (1)  will  be  rewritten  as 

x  =  G(x),  (2) 

where  G(x)  =  x  —  F(x)H(x)  for  any  H  such  that,  if  s  is  a  solution  of  (1),  then 
E r(s)  is  finite  and  nonzero  (Traub,  1964).  If  s  is  a  solution  of  (1),  then  s  is  also  a 
solution  of  (2)  and  we  have  converted  our  problem  to  determining  a  fixed  point  s 
of  Equation  (2).  A  sequence  {xn}  is  produced  from  Equation  (2)  by  the  iterative 
scheme 

xn+l  =  G(xn),  n  =  0,1,...,  (3) 

where  limn_0o  xn  =  s  for  a  convergent  sequence  {xn}.  The  function  G  is  referred 
to  as  an  iteration  function.  An  iteration  is,  therefore,  defined  as  computing  xn+1 
by  Equation  (3)  for  some  non-negative  integer  n.  An  acceleration  technique  will 
produce  a  new  sequence  {yn}  which  also  converges  to  s,  but  faster  than  the  original 
sequence  {xn}. 

Sequences  are  generated  several  different  ways.  For  scalars,  there  exist  the 
well-studied  sequences  produced  by  the  partial  sums  of  a  series.  Perhaps  the  best 
known  method  for  producing  vector  sequences  is  numerically  solving  the  system  of 
linear  equations 

x  =  Ax  +  b,  (4) 

where  A  is  an  m  x  m  matrix,  x  and  b  axe  m-dimensional  column  vectors  with  b 
constant.  Hence,  the  basic  iteration  equation  becomes 

xn+i  =  Axn  +  b=  G{xn).  (5) 

Vector  sequences  can  also  be  generated  by  nonlinear  problems:  integral  equations 
and  ordinary  differential  equations.  Examples  and  test  problems  of  varied  form  will 
be  presented  throughout  +he  thesis. 
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Equations  (4)  and  (5)  will  be  used  extensively  in  later  chapters  in  the  devel¬ 
opment  of  the  acceleration  methods.  However,  in  practice,  problems  of  this  type 
would  normally  be  solved  by  other  methods  that  will  not  be  discussed  in  this  paper. 
The  small  linear  problems,  where  the  value  of  m  in  Equation  (4)  is  small,  would 
be  solved  by  a  direct  method  called  Gaussian  elimination.  There  also  exist  efficient 
special  methods  that  will  solve  sparse  linear  problems  of  large  dimension.  However, 
linear  problems  are  quite  useful  for  designing  extrapolation  models  and  testing  the 
algorithms.  The  small  nonlinear  problems  would  be  solved,  in  practice,  by  Newton’s 
method  or  an  optimization  method,  for  example,  a  quasi-Newton  method.  These 
problems  are  most  useful  as  test  problems  where  the  limit  can  only  be  estimated 
and  not  found  exactly.  The  domain  of  problems  that  will  be  most  practical  for  the 
methods  presented  in  this  study  is  medium-to-large- sized  nonlinear  problems.  Fox 
(1965),  Ortega  and  Rheinboldt  (1970),  Varga  (1962),  and  Young  (1971)  provide 
further  background  material  on  these  other  methods. 


CHAPTER  II 


MOTIVATIONAL  EXAMPLES 
Using  Acceleration  Techniques 

Before  developing  the  acceleration  techniques,  this  chapter  will  be  used  to  help 
motivate  interest.  The  motivation  will  come  in  two  parts.  First,  it  will  be  shown 
how  an  acceleration  technique  can  be  used  to  find  the  limit  of  a  sequence.  It  will 
then  be  shown  how  the  eigenvalues  of  the  matrix  A  of  the  iteration  Equation  (5) 
effect  the  convergence  of  the  vector  sequence  {xn}. 

To  start  with,  consider  the  scalar  problem 

F(x)  =  x-  e~0  Sm  =  0,  (6) 

with  a  solution  of  0.7034674  for  seven  place  accuracy.  Figure  1  (page  6)  is  the 
graph  of  F(x)  on  the  interval  [0, 1].  The  solution  is  found  by  solving  the  fixed  point 
solution  of  (2),  where 

x  =  e-0'5*  =  G{x). 

This  will  give  an  iteration  equation  of 

*n+1  =  e-0'5*"  =  G(xn).  (7) 

Figure  2  (page  6)  shows  the  graph  of  G(x)  and  the  equation  y  =  x  on  the  inter¬ 
val  [0,1],  In  addition,  the  figure  illustrates  graphically  the  convergence  of  the  fixed 
point  problem  (Conte  and  de  Boor,  1980)  with  an  initial  'value  of  x0  =  0.5.  If  a  solu¬ 
tion  of  Equation  (6)  exists,  then  it  will  be  the  intersection  of  y  =  x  and  y  =  G(x).  To 
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Figure  2.  Graphs  of  y  =  G(x )  =  e  °  5*  and  y  =  x 


find  this  point  graphically,  it  is  known  that  the  point  (xn_i,xn)  must  lie  on  the 
graph  of  G(x).  The  next  point,  (xn,xn+i),  can  be  found  by  drawing  a  line  through 
(xn_i,xn)  parallel  to  the  x  axis.  This  line  intersects  the  graph  of  y  =  x  at  (xn,xn). 
Now  draw  through  this  point  the  line  parallel  to  the  y  axis.  The  intersection  of  this 
line  and  the  graph  y  =  G(x )  will  be  the  point  (xn,(?(xn)),  or  (xn,i„+i).  Continuing 
this  procedure  will  eventually  lead  us  to  the  desired  solution. 

It  is  a  fact,  however,  that  not  all  fixed  point  iterations  converge.  It  often 
occurs  that  a  mathematical  problem  has  a  unique  and  reasonable  solution,  but 
when  a  numerical  algorithm  is  devised  to  solve  the  problem,  like  the  fixed  point 
iteration  scheme,  the  resulting  sequence  of  approximations  diverges.  However,  there 
are  certain  criteria  that  will  insure  convergence.  It  is  not  the  intent  of  this  paper 
to  discuss  such  criteria,  but  one  can  find  good  discussion  for  the  linear  case  in 
texts  written  by  Henrici  (1964)  or  by  Conte  and  de  Boor  (1980).  We  will  consider 
divergent  sequences  after  determining  the  solution  of  Equation  (6). 

Using  (7)  and  x0  —  0.5,  a  sequence  {xn}  is  derived  whose  limit  is  the  solution. 
As  shown  in  Table  1  (page  8),  it  requires  16  iterations  before  the  sequence  converges 
to  the  correct  value  with  seven  decimal  place  accuracy.  The  correct  digits  for  each 
iterate  are  underlined.  In  addition,  Table  1  gives  the  differences  between  consecutive 
iterates  and  also  the  ratios  of  consecutive  differences,  denoted  by  un  =  xn+i  —  xn 
and  rn  =  un/un_!,  respectively.  This  information  will  help  develop  an  acceleration 
technique  to  apply  to  Equation  (7)  for  comparison  (Atkinson,  1972). 

A  closer  look  at  the  first  few  ratios  of  Table  1  shows  that  they  seem  to  be  con¬ 
verging  to  a  value  of  approximately  —0.3517.  The  ratios  for  n  >  9  no  longer  converge 
due  to  rounding  errors  in  determining  xn  and  un  to  seven  places.  Assuming  that  the 
ratio  is  approximately  constant  after  the  fifth  iteration,  estimates  for  u5  through 
un,  denoted  by  un,  n  =  5, . . . ,  11,  can  be  made  by  Uj+1  =  u,(f),  i  =  4, . . . ,  10,  where 
f  =  r4  =  —0.3512072  and  u4  =  u4.  These  estimates  are  given  in  Table  2  (page  9) 


along  with  the  actual  values  from  Table  1. 


TABLE  1 

ITERATED  VALUES  FOR  EQUATION  (7) 


n 

u  =  zn+1  -  xn 

rn  = 

0 

0.5000000 

0.2788008 

1 

0.7788008 

-0.1013378 

-0.3634773 

2 

0.6774630 

0.0352108 

-0.3474596 

3 

0.7126738 

-0.0124371 

-0.3532183 

4 

0.7002367 

0.0043680 

-0.3512072 

5 

0.7(146347 

-0.0015372 

-0.3519230 

6 

0.7030675 

0.0005406 

-0.3516783 

7 

0.7036081 

-0.0001902 

-0.3518312 

8 

0.7034179 

0.0000669 

-0.3517350 

9 

0.7034848 

-0.0000235 

-0.3512705 

10 

0.7034613 

0.0000083 

-0.3531914 

11 

0.7Q34696 

-0.0000029 

-0.3493975 

12 

0.7Q34667 

0.0000010 

-0.3448275 

13 

0.7034677 

-0.0000004 

-0.4000000 

14 

0.7034673 

0.0000002 

-0.5000000 

15 

0.IQMSI5 

-0.0000001 

-0.5000000 

16 

0.7034674 

0.0000000 

0.0000000 

17 

0.7034674 
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TABLE  2 

ESTIMATED  AND  ACTUAL  VALUES  FOR  CONSECUTIVE 
DIFFERENCES  FOR  EQUATION  (7) 


n 

un  =  un. 

-i(f) 

«n(r) 

5 

(0.0043680)f  = 

-0.0015340 

-0.0015372 

6 

(— 0.0015340)f  : 

=  0.0005388 

0.0005406 

7 

(0.0005388)f  = 

-0.0001892 

-0.0001902 

8 

(— 0.0001892)f 

=  0.0000664 

0.0000669 

9 

(0.0000664)f  = 

-0.0000233 

-0.0000235 

10 

(— 0.0000233)f 

=  0.0000082 

0.0000083 

11 

(0.G000082)f  = 

-0.0000029 

-0.0000029 

It  is  true  that 

*ii  =  *5  +  (*«  —  *5)  4-  (*7  —  *e)  4-  •  •  •  +  (zn  —  X10) 

=  x5  4-  u5  4-  4-  u7  +  us  +  u9  4-  Ui0. 

Therefore,  we  can  estimate  in  by 

x'u  «  0.7046047  -  0.001534  +  0.0005388  -  0.0001892  + 

0.0000664  -  0.0000233  4-  0.0000082  -  0.0000029 
=  0.7034687. 

Hence,  using  only  information  obtained  from  X3,  X4,  and  X5;  x'n  has  been  determined 
more  accurately  than  the  actual  iterated  xn.  Even  though  an  estimate  for  xj,,  for 
some  positive  integer  k ,  may  not  always  be  more  accurate  than  the  iterated  Xfc,  the 
estimate,  x'k,  will  usually  be  a  better  approximation  of  the  solution  than  xp,  for 
some  p,  p  <  k. 
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Generalizing  the  concept,  assume  zn_2,  xn-i>  and  xn  have  been  computed. 
Also  assume  rn  is  approximately  constant  with  f  =  rn_i  =  un_i/un_2.  Then  the 
solution  can  be  approximated  by  viewing  it  as  the  limit,  z^,,  of  the  sequence: 


xoo  ~  +  'Un  +  ^n+1  +  Un+2  +  '  *  ' 

~  Xn  +  un-i f  +  unf  +  un+if  H - 

~  xn+  Un-l f  +  Un_if2  +  Un_if3  H - 

=  xn  +  +  r2  +  f3  H - ) 

=  z„  +  ttn_i  (f/(l  -  f)) ,  |f|  <  1,  (8) 

since  the  series  in  parentheses  is  a  geometric  series.  Substituting  f  =  un-i/un- 2 
and  simplifying  give 


Xoo  ^  Xn  un— 1  ( [un-l  /^n-j]  /  [(ttn-l/^n-2 )  1]) 

=  Xn  ('Un-l)  /(tin- 1  “  tin-2) 

=  *»  - 


(z„  -  Zn_x)2 


(9) 


(*n  —  ®n-l)  —  (*n-l  ~  *n-2) 

This  formula  is  Aitken’s  A2  formula  (1936-37)  for  accelerating  a  convergent  se¬ 
quence.  More  information  will  be  discussed  concerning  Aitken’s  formula  in  Chapter 
IV.  However,  based  on  the  above  development,  one  may  assume  that  if  Aitken’s 
method  is  applied  to  Equation  (7)  after  the  nth  iteration,  a  better  estimate  of  the 
solution  can  often  be  found  with  no  additional  iterations. 

Aitken’s  formula  can  be  used  as  a  sequence  generator  to  derive  a  new  sequence 
{yn}.  The  new  sequence  is  generated  by  Formula  (9)  rewritten  as 

(®n+2  ~  <cn+l)2 


Vn+2  —  *n+2 


xn+2  2zn+J  4-  Xn 


,  71  =  0,1,... 


(10) 


In  fact,  (10)  can  be  used  to  produce  even  a  third  sequence  {zn}  from  {y„}.  Apply¬ 
ing  this  technique  to  the  sequence  derived  from  Equation  (7),  the  resulting  three 
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sequences  are  shown  in  Table  3.  One  can  see  that  the  solution,  to  seven  place  ac¬ 
curacy,  is  found  after  only  five  iterations.  Generating  sequences  in  this  fashion  is 
called  a  static  model  and  Table  3  is  referred  to  as  a  static  display. 


TABLE  3 

AITKEN’S  COMPUTED  VALUES  FOR  EQUATION  (7) 


n 

Xn 

Vn 

0 

0.5000000 

1 

0.7788008 

2 

0.6774630 

0.1044777 

3 

0.7126738 

0.7035942 

4 

0.7002367 

0.7034830 

0.7034669 

5 

0.7046047 

0.7034693 

0.7034674 

Equation  (10)  is  written  differently  than  the  way  it  is  given  in  many  texts. 
The  difference  is  that  the  new  term  in  {y„}  is  denoted  by  the  subscript  n  +  2  instead 
of  its  usual  subscript  n.  The  reason  for  this  change  is  that  if  we  want  to  compare 
terms  of  the  two  sequences  for  error,  that  is,  their  closeness  to  the  exact  answer, 
then  to  compare  the  nth  term  of  {yn},  call  it  y  for  the  moment,  to  the  nth  term  of 
{in},  call  it  x,  is  unfair.  This  y  cannot  be  computed  until  xn+2  is  available.  If  the 
original  sequence  converges,  not  only  should  y  have  a  smaller  error  than  x,  but  so 
should  the  next  two  terms  following  x  in  the  original  sequence.  Therefore,  what  is 
valuable  is  to  compare  y  with  xn+2  in  error.  Hence,  the  nth  term  of  {yn}  is  referred 
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to  as  yn+2  so  that  when  comparisons  are  made,  the  elements  being  compared  will 
have  the  same  subscripts. 

Now  consider  the  geometric  series 

1  +  2  +  4  +  . . . 

Then  the  sequence  of  partial  sums  is 

1,3,7,...,  or  xn  =  — 1  +  2n,  n=l,2,...  (11) 

If  Aitken’s  method  is  applied  to  Sequence  (11),  the  resulting  sequence  is 

-1,-1, -1,... 

Hence,  the  acceleration  technique  determined  the  value  —1  as  the  “limit”  of  a 
divergent  sequence.  Shanks  (1955)  says  the  divergent  sequence  is  “diverging  from” 
—1  and  calls  the  value  —1  the  “antilimit”  of  (11). 

As  stated  earlier,  iterated  sequences  of  some  mathematical  problems  do  not 
converge.  However,  if  the  sequence  has  an  antilimit,  the  antilimit  is  usually  unique, 
equal  to  the  solution  of  the  problem,  and  can  usually  be  found  by  applying  an  ac¬ 
celeration  method  to  the  original  divergent  sequence  (Sidi,  Ford,  and  Smith,  1986). 
There  are  cases  known  where  this  is  not  true  (Shanks,  1955);  therefore,  one  must  be 
careful  to  ensure  that  the  computed  antilimit  is,  in  fact,  the  solution.  An  example 
where  Aitken’s  method  gives  erroneous  results  will  be  discussed  in  Chapter  IV. 

The  early  leader  in  the  use  of  divergent  sequences  to  derive  correct  answers 
was  Euler  (1707-1783).  He  maintained  that  if  a  function  /  gave  rise  to  a  series,  then 
the  “sum”  of  the  series  should  be  f(x)  for  any  x,  even  when  the  series  diverged. 
Even  though  his  definition  of  the  word  “sum”  extended  the  normal  definition  of 
a  sum  of  a  series,  he  felt  quite  comfortable  with  it  since  “the  new  definition  ... 
coincides  with  the  ordinary  meaning  when  a  series  converges...”  (Bromwich,  1926, 
p.  322). 
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One  may  look  at  antilimits  as  the  assigning  of  a  number  to  a  divergent  se¬ 
quence.  This  has  been  around  for  quite  some  time  in  mathematics  in  the  form  of 
summability  methods.  There  exists  several  methods  of  summability.  Excellent  ma¬ 
terial  on  the  subject  matter  may  be  found  in  texts  written  by  Hardy  (1949),  Lubkin 
(1952),  Moore  (1938),  and  Zygmund  (1959). 

Illustrations  of  Convergence  in  Two  Dimensions 

Let  us  now  look  at  the  second  motivational  factor  of  this  chapter.  Consider 
Equation  (4)  with  dimension  m  =  2.  Define  A  and  b  by 

1.0  0.1  _  1.2 

A  =  and  b  =  .  (12) 

-0.5  0.4  -2.0 

Using  Equation  (5)  and  the  initial  vector  (1,1),  a  convergent  sequence  {xn}  of  two- 
dimensional  vectors  is  generated  that  converges  to  the  solution  s  =  (10.4,  — 12.0)r. 
Table  4  shows  the  first  six  terms  of  the  sequence.  The  path  of  convergence  of {x„} 
to  j*is  shown  in  Figure  3  (page  14).  The  graph  can  be  considered  as  a  trajectory  of 
a  moving  particle  originating  at  x0  and  terminating  at  the  solution. 


TABLE  4 

FIRST  SIX  TERMS  OF  THE  GENERATED 
SEQUENCE  OF  PROBLEM  (12) 


n 

xn 

n 

Xn 

0 

(1.000000,  1.000000) 

3 

(4.091000,  -5.241000) 

1 

(2.300000,  -2.100000) 

4 

(4.766900,  -6.141900) 

2 

(3.290000,  -3.990000) 

5 

(5.352710,  -6.840210) 
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Figure  3.  Trajectory  for  Eigenvalues 
0,9  and  0.5 


The  matrix  A  has  eigenvalues  q  =  0.9  and  q\  =  0.5  with  associated  eigenvectors 
v  =  (1,  —  l)r  and  vi  =  (1,  —  5)T,  respectively,  which  are  also  shown  in  Figure  3.  The 
trajectory  of  {*„}  shows  that  as  xn  converges  to  /,  the  convergence  is  asymptotic 
along  the  vector  v.  Since  v  and  are  linearly  independent,  the  vector  x0  —  /  can 
be  written  as  a  linear  combination  of  these  vectors.  Hence, 

x0-s  =  (-9.4, 13.0)r  =  -8.5(1,  -1)T  -  0.9(1,  -5)t. 

=  — 8.5u  —  0.9i/i. 

Using  the  fact  that  v  and  Vi  are  eigenvectors,  it  follows  that 

*i  -  /  =  A{x o  -  s)  =  -8.5(0.9)v  —  0.9(0.5)Fi 
=  (2.3,  -2.1)r,  and 

x2-s  =  A(x j  -  s)  =  — 8.5(0. 9)2u  —  0.9(0.5)2u1 
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=  (3.29,  -3.99)r. 

Continuing,  it  can  be  seen  that  the  sequence  {xn}  can  be  generated  by  the  equation 

xn  =  s  +  aqnv  +  ai5"ui, 

where  a  =  —8.5  and  a\  =  —0.9.  Assume  that  n  =  100.  Then 

x100  =  i'+a(0.9)lootr+a1(0.5)loV1 

w  5*+  (2.66  x  10~5)au  +  (7.89  x  lO-31)^^. 

Since  the  last  term  is  approximately  zero,  z10o  is  primarily  the  sum  of  s  and  a 
multiple  of  the  eigenvector  v.  So,  as  n  increases,  the  convergence  of  the  sequence  is 
controlled  by  the  dominant  eigenvadue,  0.9,  resulting  in  the  asymptotic  convergence 
along  v.  In  addition,  for  n  greater  than  200,  the  coefficient  of  v  is  less  than  1  x  10-8 
which  implies  {afn}  has  converged  to  s  with  seven  place  accuracy. 

The  problem  can  be  generalized  for  dimension  m  where  q1,q2, . . .  ,qm  are  the 
eigenvadues  of  A  with  corresponding  eigenvectors  v\,  v2, . . . ,  vm;  and  with  the  as¬ 
sumption  that  unity  is  not  an  eigenvalue  of  A  so  that  the  problem  has  the  unique 
solution  i*.  Then  for  some  scalars  a;,  the  sequence  {xn}  can  be  generated  by 

m 

*n  =  s- +  a*  ,  n  =  0,l,...  (13) 

«=i 

If  qi  is  the  dominant  eigenvadue  and  ai  ^  0,  which  will  usually  be  true  for  the  given 
£o,  then  the  limit  of  {zn}  is  5*  provided  |gi|  <  1.  If  |?i|  >  1,  then  {®n}  is  a  divergent 
sequence  and  i*  is  its  antilimit.  Therefore,  assuming  that  {zn}  converges,  it  is  the 
dominant  eigenvalue,  call  it  q,  that  is  of  importance.  The  smaller  the  modulus  of  q, 
the  faster  the  convergence.  However,  if  one  or  more  eigenvadues  have  modulus  close 
to  unity,  then  the  convergence  will  be  slow. 

Figures  4  through  6  (pages  16  and  17)  show  other  trajectories  of  the  sequence 
{x„}  generated  by  Equation  (5)  with  dimension  two.  The  matrix  A  and  associated 
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eigenvalues  for  each  figure  are 

Figure  4  Figure  5  Figure  6 


-1.0 

-0.2 

A  = 

1.0 

-0.2 

A  = 

1.0 

-1.0 

0.75 

0.6 

0.75 

-0.6 

1.0 

-0.2 

*  =  -0.9, 0.5  0.9, -0.5  $  =  0.4  ±  0.8* 

In  addition,  the  eigenvectors  are  graphed  and  labeled  as  v,  for  the  eigenvector 
corresponding  to  the  dominant  eigenvalue  q  and  referred  to  as  the  major  eigen¬ 
vector,  and  vu  the  eigenvector  corresponding  to  qi  and  referred  to  as  the  mi¬ 
nor  eigenvector.  Since  the  eigenvalues  of  the  matrix  A  in  Figure  6  are  complex 
conjugates,  the  eigenvectors  are  also  complex  and,  hence,  are  not  graphed.  Fig¬ 
ures  4  and  5  support  the  fact  that  the  convergence  of  {:?„}  is  indeed  asymp¬ 
totic  along  the  eigenvector  v.  Figure  6  suggests  that  the  trajectory  of  {afn}  for 
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Figure  5.  Trajectory  for  Eigenvalues 
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Figure  6.  Trajectory  for  Eigenvalues 
0.4  ±  0.8i 
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problems  with  complex  eigenvalues  is  some  form  of  a  spiral.  The  trajectory  of 
another  problem  with  complex  conjugate  eigenvalues  is  shown  in  Figure  7. a.  Re¬ 
turning  to  Figure  3,  it  shows  a  smooth  monotonic  convergence.  This  is  not  always 
the  case  for  two  positive  eigenvalues.  However,  as  n  continues  to  increase  the  conver¬ 
gence  will  eventually  become  monotonic  and  resemble  Figure  3.  Figure  7.b  shows 
the  trajectory  of  another  problem  with  positive  eigenvalues  that  initially  begins 
to  “converge”  along  the  minor  eigenvector  but  eventually  converges  monotonically 
along  the  major  eigenvector.  Figure  7.c  shows  an  example  of  a  trajectory  for  a 
problem  with  two  negative  eigenvalues. 

Comparing  Figures  4  and  5,  we  see  that  both  trajectories  zig-zag  across  the 
eigenvector  associated  with  the  negative  eigenvalue.  Because  this  is  the  major 
eigenvector  in  Figure  4,  the  zig-zag  motion  dampens  out  as  the  sequence  approaches 
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Figure  7.  Other  Trajectory  Examples 
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s,  and  the  convergence  becomes  almost  linear  along  the  major  eigenvector.  Figure 
7.d  shows  the  trajectory  of  a  problem  where  the  case  is  reversed  and  the  graph 
starts  in  an  almost  linear  approach  and  eventually  ends  up  in  the  zig-zag  motion. 
Why  do  the  trajectories  of  these  two  problems  differ? 

From  Equation  (13),  we  have 

m 

xn  —  3  =  Oj  Vj  q ”,  n  =  0, 1, . . . 

t=i 

Since  the  left-hand  side  of  the  equation  is  the  difference  between  the  nth  term  of 
the  sequence  {xn}  and  the  solution  s,  it  will  be  referred  to  as  the  nth  error  vector, 
denoted  by  en.  For  m  =  2,  the  error  vector  can  be  rewritten  as 

en  =  axvxq ”  +  avqn, 

where  q  is  the  dominant  eigenvalue  and  v  is  the  eigenvector  associated  with  q.  If 
n  =  0,  then  the  initial  error  vector  is 

e*o  =  a  iFi  +  av. 

If  a  is  sufficiently  smaller  than  ax,  then  the  dominant  eigenvector  for  the  first  few 
sequence  elements  will  be  t?i.  Therefore,  the  trajectory  will  start  almost  parallel 
to  vx.  However,  as  n  increase,  the  higher  powers  will  result  in  a  dampening  of  qx 
and  the  dominant  eigenvalue  will  cause  a  convergence  in  the  direction  of  the  major 
eigenvector.  Figure  4  shows  the  resulting  trajectory  for  q  negative  and  qx  positive. 
If  the  signs  of  the  eigenvalues  are  reversed,  then  the  trajectory  starts  almost  linearly 
and  terminates  in  a  zig-zag  fashion  as  shown  in  Figure  7.d. 

Relaxation  Factor 

Up  to  this  point,  successive  elements  of  the  generated  sequence  {£„}  have  been 
found  by  using  the  iteration  Equation  (3).  However,  there  is  a  variation  of  (3)  that 
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may  be  helpful  in  accelerating  the  convergence.  This  method  uses  a  constant  called 
a  relaxation  factor  to  adjust  the  distance  the  iteration  moves  from  the  previous 
sequence  element.  The  new  iteration  equation  is 

yn+ 1  =  Gw{yn)  =  yn  +  w(G{yn)  -  yn),  (14) 

where  n  =  0, 1, . . . ,  yo  =  xo  °f  Equation  (5),  and  w  >  0  is  the  relaxation  factor.  For 
0  <  w  <  1,  (14)  is  called  vector  under-relaxation,  and  for  w  >  1,  (14)  is  referred  to 
as  vector  over-relaxation.  With  no  restrictions  on  w,  (14)  is  a  parametric  equation 
for  the  line  containing  the  points  yn  and  G(yn).  For  w  =  0  and  1,  yn+i  equals  yn 
and  G(yn),  respectively.  If  0  <  w  <  1,  yn+1  is  located  on  the  line  between  y  and 
(?(^n).  If  w  >  1,  yn+i  is  still  on  the  line;  however,  yn+ 1  is  “beyond”  G(yn),  as  viewed 
form  yn. 

Returning  to  (14)  and  using  Equation  (5), 

yn+ 1  =  yn  +  w([Ayn  +  b]-yn) 

=  w(Ayn  +  b)  +  (l-  w)yn 

=  (wA  -I-  [1  -  w]I)yn  +  wb 

=  Awyn  -f  wb,  (15) 

where  Aw  =  wA  +  (1  —  w)I.  Equations  (14)  and  (15)  axe  equivalent  iteration 
equations. 

Let  qi  and  pi ,  i  =  l,...,m,  be  the  eigenvalues  of  the  matrices  A  and  Aw, 
respectively.  Then 

0  =  det(Au,  —  pi)  =  det(iu.4  +  [1  —  w]I  —  pi) 

=  det(u>i4  +  [1  —  w  —  p]I) 

—  u;mdet(>l  —  [l  —  (1  —  p)/w\I). 

Since  w  ^  0,  the  eigenvalues  of  A  must  be  q  =  l  —  (l—p)/w.  Hence,  p  =  l+w(q  —  1). 
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Therefore,  we  have  a  parametric  equation  for  a  line  through  the  point  z  =  1  =  (1, 0) 
and  q  in  the  complex  plane.  Since 

\P-  1|  =  I1  +*»(?  -  1)  ~  1|  =  \w(q  -  1)|  . 

the  distance  between  the  eigenvalues  of  Aw  and  the  point  z  is  w  times  the  distance 
between  the  eigenvalues  of  the  matrix  A  and  z. 

As  a  result,  Equation  (13)  can  be  rewritten  with  a  new  set  of  eigenvalues. 
By  carefully  choosing  w,  one  may  be  able  to  convert  a  divergent  problem  into  a 
convergent  one.  For  example,  if  we  have  a  two-dimensional  problem  with  eigenval¬ 
ues  0.5  and  —1.5,  then  iteration  Equation  (5)  will  generate  a  divergent  sequence. 
However,  the  eigenvalues  of  iteration  Equations  (14)  and  (15)  with  a  relaxation 
factor  of  w  =  0.5  are  0.75  and  —0.25.  Hence,  by  using  Equation  (14),  one  “may” 
produce  a  convergent  sequence  even  if  some  of  the  eigenvalues  of  A  have  moduli 
greater  than  unity.  The  word  may  is  used  because  there  are  cases  where  no  value  of 
w  will  convert  a  divergent  sequence  into  a  convergent  one;  for  example,  a  problem 
with  an  eigenvalue  of  2.0.  For  any  value  of  w  ^  0,  the  modulus  of  the  converted 
eigenvalue  will  always  be  greater  than  unity.  Figure  8  (page  22)  shows  the  region 
in  the  complex  plane  of  eigenvalues  of  Equation  (4)  for  which  Equation  (14)  with 
w  —  0.5  generates  a  convergent  sequence.  The  region  can  be  described  as 

{p  :  |p  —  d|  <  4,  where  o  =  —  1  +  0i}. 

Choosing  w  >  1  can  also  be  useful  on  some  occasions.  Consider  the  two- 
dimensional  problem  where  the  eigenvalues  are  0.6  and  0.9.  Even  though  the  itera¬ 
tion  equation  will  produce  a  convergent  sequence,  because  the  dominant  eigenvalue 
is  close  to  unity,  the  convergence  will  be  slow.  A  relaxation  factor  of  w  =  2  will 
transform  the  problem  into  an  equivalent  problem  with  eigenvalues  0.2  and  0.8. 
Hence,  the  transformation  has  a  smaller  dominant  eigenvalue,  resulting  in  faster 
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Figure  8:  Eigenvalue  Region  of  A  Resulting 
in  Convergent  Sequences  for 
Aw,  w  =  0.5 


convergence. 

One  must  show  care  in  choosing  w.  There  are  problems  where  a  relaxation 
factor  less  than  unity  can  cause  a  convergent  sequence  to  converge  slower  or  a  relax¬ 
ation  factor  greater  than  unity  can  cause  a  convergent  sequence  to  diverge.  Consider 
the  sequence  generated  by  (5)  with  eigenvalues  0.8  and  —0.2.  Using  Equation  (14) 
and  w  =  0.5,  the  transformation  problem  has  eigenvalues  0.9  and  0.4.  The  trans¬ 
formed  problem  has  a  dominant  eigenvalue  closer  to  unity;  hence,  the  convergence 
is  slower.  If  w  =  2.0,  the  transformation  problem  has  eigenvalues  0.6  and  —1.4, 
resulting  in  a  divergent  problem.  Unless  otherwise  stated,  the  results  given  in  this 
study  will  be  for  w  —  1. 


CHAPTER  III 


GENERALIZED  INVERSES 

Before  we  continue  with  acceleration  methods,  there  is  an  area  that  needs 
extended  coverage  beyond  that  which  is  available  in  a  normal  course  of  matrix 
theory  or  linear  algebra.  This  area  is  the  theory  and  applications  of  what  is  called 
a  generalized  inverse  of  a  matrix. 

Let  A  be  a  square  m  x  m  matrix  with  rank  R(A)  =  m.  Then  we  know  that 
there  exists  a  unique  matrix  B,  called  the  inverse  of  A ,  such  that  AB  =  BA  =  Im, 
where  Im  is  the  identity  matrix  of  order  m.  The  inverse  of  A  is  normally  denoted 
by  A-1 .  Hence,  given  the  square  mxm  matrix  A  and  the  m-dimensional  vector  y , 
then  the  solution  of  a  set  of  consistent  linear  equations 

y  =  Ax  (16) 

is  (A_1)y  =  (A_1)(Ax)  =  Imx  =  x.  If  A  has  an  inverse,  then  A  is  said  to  be 
nonsingular;  otherwise,  A  is  singular.  If  A  is  a  rectangular  matrix,  then  no  such 
matrix  B  exists  as  the  inverse  of  A  and,  thus,  a  simple  expression  of  a  solution  of 
(16)  in  terms  of  A  is  more  difficult. 

Moore  (1920)  extended  the  normal  concept  of  inverses  to  singular  and  rectan¬ 
gular  matrices.  However,  the  theoretical  properties  of  these  matrices  were  not  fully 
investigated  until  1955,  when  Penrose  defined  a  uniquely  determined  inverse  matrix 
for  any  matrix  A  which  he  called  the  generalized  inverse.  Moore’s  and  Penrose’s 
inverses  are  equivalent  when  the  inner  product  of  the  two  m-dimensional  vectors 
x  =  (*!,.. .  ,xm)T,y  =  (yi,..  .,ym)T  is  defined  by 


23 


24 


(*>y)  =  (»*)*  =  X)y;  x»’ 

i=i 

where  y  *  indicates  the  complex  conjugate  transpose  of  y  and  j/j  is  the  complex 
conjugate  of  y,.  Penrose’s  definition  of  a  generalized  inverse  is  as  follows: 


Definition:  For  any  matrix  A,  square  or  rectangular,  real  or  complex,  there  exists  a 
unique  matrix  G  satisfying  the  conditions 


(1)  AG  A  =  A  (2)  GAG  =  G 

(3)  (AG)'  =  AG  (4)  {GA)*  =  GA. 

G  is  called  the  Moore-Penrose  generalized  inverse  of  A. 


(17) 


With  the  use  of  generalized  inverses,  we  can  extend  the  concept  of  solving  (16) 
where  A  is  a  m  x  k  matrix  of  rank  r,  r  —  k  <  m,  y  is  a  m-dimensional  vector,  and 
x  is  a  A:-dimensional  vector.  If  r  =  k  =  m,  A  is  nonsingular  and  the  problem  is 
as  before  with  a  solution  of  x  =  (A-1)y.  If  m  ^  k,  then  x  =  Gy,  where  G  is  the 
unique  solution  of  Equations  (17).  However,  as  Penrose  pointed  out,  a  solution  of 
(16)  does  no  require  a  matrix  G  which  satisfies  all  the  conditions  of  (17).  One  can 
find  a  solution  of  (16)  which  satisfies  only  condition  (1)  of  (17).  Given  y  =  Ax,  then 
AGy  =  AG  Ax  =  Ax.  Therefore,  x  =  Gy. 

Some  authors  refer  to  generalized  inverses  by  other  names.  Greville  (1959)  and 
Rohde  (1964)  prefer  the  use  of  the  name  “pseudo-inverses.”  Rao  (1965)  referred  to 
them  frequently  as  the  “Moore-Penrose  inverses.”  Albert  (1972)  combines  the  names 
together  and  calls  them  the  “Moore-Penrose  pseudo-inverse.”  In  addition,  different 
names  are  given  to  matrices  which  satisfy  one  or  more  of  the  conditions  of  (17). 
In  this  study,  generalized  inverses  will  refer  to  those  matrices  that  satisfy  at  least 
condition  (1)  of  (17). 
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Generalized  inverses  of  matrices  satisfying  condition  (1)  of  (17)  are  not  unique 
(Pringle  and  Rayner,  1971).  Let  A  be  a  m  x  k  matrix.  If  R(A)  =  m,  then  a 
generalized  inverse  of  A  is 

G  =  {V  A*)(AV  A*)~l , 

where  V  is  an  arbitrary  matrix  such  that  R{AV A*)  =  i2(A).  G  is  called  a  right 
inverse  of  A  in  this  case.  If  R(A)  =  k,  then  a  generalized  inverse  of  A  is 

G  =  ((A*)FA)-1(A*)F,  (18) 

where  V  is  an  arbitrary  matrix  such  that  R{{A*)V A)  —  R{A).  G  is  called  a  left 
inverse  of  A  for  this  case. 

Returning  to  the  problem  (16),  when  A  is  a  m  x  k  matrix  of  rank  r  =  k  <  m, 
we  can  find  a  solution  to  the  problem  by  using  generalized  inverses.  Hence, 

x  =  Irx  =  GAx  =  Gy 

where  G  is  of  the  form  (18).  A  simple  choice  for  V  is  V  =  Ir  such  that 

G  =  (A*A)~1At.  (19) 

Henceforth,  the  definition  of  the  generalized  inverse  of  the  m  x  k  matrix  A,  k  <  m, 
is  the  k  x  m  matrix  G  as  defined  in  (19)  with  the  notation  G  =  A+. 


CHAPTER  IV 


AITKEN’S  A2  METHOD  FOR 
SCALARS  AND  VECTORS 

Theory  for  Scalars 

Aitken’s  A2  method  was  introduced  in  Chapter  II  to  help  solve  Equation  (6). 
Consider  Equation  (8)  applied  to  a  geometric  series  where  {xn}  is  the  sequence  of 
partied  sums.  Then  the  ratio,  r,  is  constant  and,  hence,  (8)  can  be  written  as  the 
sum 

s  =  xn  +  un_i  ( ,  |r|  <  1. 

Substituting  u„_i  =  run_2  =  rn-1u0  and  u0  =  =  rx0  gives 

3  =  +  (j^) 

=  *- +  (rr?) (20) 

Therefore,  for  |r|  <  1,  xn  —  s  =  crn,  where 

—  XqT 

°  ~  1  -  r* 

So,  there  exist  constants  s  and  r  such  that  as  n  increases  by  one,  the  distance  from 
xn  to  3  is  multiplied  by  r.  It  is  easy  to  see  that  if  |r|  <  1,  then  rn  goes  to  zero  as 
n  increases,  which  implies  that  s  is  the  limit  of  the  sequence  {®n}-  If  r  =  —  1  or 
[r|  >  1,  then  s  is  the  antilimit  of  {in}-  This  is  a  special  case  of  Equation  (13). 

If  we  assume  that  there  exists  a  constant  p  ^  s  such  that  xn  —  p  =  drn,  for 
some  constant  d ,  then 
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drn  =  xn  —  p  =  s  +  crn  —  p. 


Hence,  s  —  p  =  (d  —  c)rn  =  brn.  Since  s  and  p  are  constants,  brn  must  remain 
constant  as  n  increases;  hence,  6  =  0.  So  s  —  p  =  0,  which  implies  s  =  p.  Therefore, 
s  is  unique. 

Written  below  is  (10)  in  a  slightly  different  form  and  two  variations: 

(*n+l)2 


Vn+2  —  Zn+2  ~ 


K) 

xnxn+  2  ~  ®n+l 


(*n+  2  —  *n+l)  —  (*n+l  ~  xn) 

K)2 


=  Xn  — 


M  ’ 


(21) 

(22) 

(23) 


where  un  =  in+1  —  xn  and  vn  =  (a:n+2  —  xn+i)  —  (®n+i  —  xn)  axe  defined  as  the 
forward  difference  operators.  Formula  (21)  is  the  most  desirable  one  for  a  convergent 
sequence  using  floating  point  arithmetic  since  zn+2  is  a  better  estimate  of  the  limit 
than  in  and  the  computed  error,  zn+2  —  s,  is  smaller  than  the  computed  error, 
xn  —  s.  On  the  other  hand,  if  s  is  the  antilimit,  then  (23.)  should  be  the  choice. 
Aitken  (1936-37)  used  (22),  but  apparently  never  used  an  automated  computer  and 
could  monitor  rounding  errors  “visually”  at  each  step. 

If  the  sequence  {zn}  is  such  that  the  ratios  of  consecutive  errors  converge  to  a 
nonzero  constant,  independent  of  n,  then  (20)  still  holds,  where  r  is  the  limit  of  the 
ratios.  Therefore,  Equations  (21)  through  (23)  will  also  hold  as  an  estimate  for  s 
and  we  can  apply  Aitken’s  method  to  approximate  the  limit  of  the  sequence.  This 
amounts  to  having  what  is  called  linear  convergence  if  |r|  <  1. 


DEFINITION:  Given  the  convergent  sequence  {xn}  and  its  limit  s.  If 

lim  ^  =  C, 
n— '°°  en 
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where  en  =  xn  —  s  is  the  error  of  the  nth  element  and  \C\  <  1,  then  the  sequence 
{xn}  converges  linearly  to  s. 


The  question  now  arises  as  to  what  happens  if  the  ratios  of  consecutive  errors 
approach  a  C  outside  of  the  open  interval  (  —  1,1).  It  can  be  shown  that  if  C  =  1, 
then  the  denominators  of  Equations  (21)  through  (23)  go  to  zero.  Hence,  the  A2 
method  will  not  work  well  if  the  ratios  converge  to  1.  If  C  is  —1  or  outside  the 
closed  interval  [—1,1],  we  do  not  have  convergence,  but  the  method  can  still  be 
used  in  hopes  of  finding  the  antilimit  of  the  sequence. 

It  was  shown  in  Chapter  II  how  Aitken’s  A2  method  can  be  used  in  the  static 
sense  to  accelerate  a  scalar  sequence.  However,  a  modification  of  this  method  is  to 
compute  the  extrapolated  value  x2  by  applying  Aitken’s  method  to  Xo,ii,  and  x2. 
Using  x'2  as  the  initial  value,  generate  two  iterates  x'3  and  x'A.  Extrapolating  again, 
we  determine  the  value  x'[  and  continue  the  pattern  until  convergence.  A  procedure 
following  this  type  of  pattern:  iterating,  extrapolating,  and  then  iterating  the  result, 
is  referred  to  as  a  repeated  method  or  a  semi-dynamic  extrapolation  model.  Figure  9 
(page  29)  shows  a  diagram  of  a  semi-dynamic  procedure  where  EXT(xn,  xn+i,  ^n^) 
implies  applying  the  extrapolation  method  to  the  values  in  parentheses. 

Vector  Theory 


For  a  sequence  of  vectors,  Aitken  (1936-37)  applied  the  extrapolation  technique 
componentwise.  In  other  words,  if  x^l\i  =  1  ,...,m,  represents  the  ith  component 
of  the  vector  x,  then  (10)  becomes 


u(,)  -x(i) 

II  n+2  ~  x  n+2 


(g(Q  V 

\X  n+2  Xn+l) 


*(0 

n+2 


x  -  2x  +  xW 


n  =  0,1,... 


The  major  problem  with  this  technique  is  that  the  computation  of  one  or  more 
component  elements  may  involve  a  denominator  of  zero,  which  obviously  results 
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Xo 

*1 

x2  EXT(x0,  xi,  x2)  ►  X2 

X3 

x\  EXT(x2,  X3,  x\)  — ►  X4  etc. 
Figure  9.  Diagram  of  Semi- Dynamic  Model 


in  an  invalid  element  in  our  new  sequence.  In  addition,  since  the  convergence  of 
the  components  of  a  vector  is  usually  related  to  one  another,  aoolying  the  tech¬ 
nique  componentwise  loses  this  relationship  and  may  result  in  some  components 
accelerating  the  “wrong”  way,  in  the  direction  opposite  to  the  majority. 

Jennings  (1971)  modified  Aitken’s  vector  method  for  the  linear  case  to  help 
prevent  the  possibility  of  this  infinite  or  even  erratic  result.  He  used  a  vector  w  to 
help  define  a  rate  of  decay  between  xn+i  and  xn+2.  His  iterative  method  is 


y  =  xn+2  +  Q(x„+2  -  xn+i),  n  =  0,l,...,  (24) 


where 


Q  = 


(w*)(xn+i  -  xn+2) 
(■UJ*)(xn+2  -  2xn+1  +  xny 


(25) 


w*  represei&s  the  complex  conjugate  transpose  of  w ,  and  y  represents  either  yn+2 
of  Equation  (10)  or  x  'n+2  depending  upon  whether  the  model  used  is  static  or 
semi-dynamic. 

In  order  to  discuss  possible  choices  of  w*  Jennings’  work,  some  definitions  are 
needed.  Because  a  m- dimensional  vector  consists  of  m  components,  it  is  convenient 
to  have  some  method  of  determining  its  size.  This  measurement  is  provided  by 
assigning  to  a  vector  a  real- valued,  nonnegative  number  known  as  a  norm.  However, 
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the  assignment  is  not  unique  since  there  exist  several  norms.  For  purposes  of  this 
study,  only  two  norms  will  be  used:  the  infinity  norm  (oo-norm)  and  the  Euclidean 
norm  (2-norm).  Given  the  vector  x  =  (xi,. . . ,  xm)T,  the  two  norms  are  respectively 
defined  as 

Hxlloo  =  max(xj|,  i  =  l,...,m;  and 

/  m  \  1/2 

11*11*  =  (£(*<  *»')2j  • 

Jennings  gave  two  choices  for  w.  The  first  choice  is  when  the  iterative  sequence 
is  governed  by  a  symmetric  matrix  A.  For  this  case,  Jennings  suggested  the  vector 
w  =  xn  —  xn+i  and  referred  to  the  formula  as  First  Difference  Modulation  (FDM). 
Given  the  scalars  a^,i  =  1, . . .  ,m,  such  that 

m 

Xo  $  =  y  ^  Of  Vi, 

»=1 

where  V{  is  the  eigenvector  corresponding  to  the  eigenvalue  qi  of  the  iterative  matrix 
A,  he  showed  that 

Q _  9»  zi 

v”E£o(i-#)*’ 

where  Z{  =  g^n(l  —  qi)2aj  >  0.  Since  all  eigenvzdues  of  A  must  have  moduli  less  than 
unity  for  convergence,  Q  cannot  have  a  zero  denominator  except  when  convergence 
is  obtained. 

Jennings’  second  choice  for  iv,  Second  Difference  Modulation  (SDM),  is  when 
the  sequence  is  governed  by  a  nonsymmetric  matrix.  For  this  case  he  chose 
w  =  xn+2  —  2x„+i  +xn.  Thus  the  denominator  is  the  square  of  the  Euclidean  norm  of 
the  second  difference  vector.  Hence,  the  denominator  cannot  result  in  zero  unless, 
once  again,  convergence  is  obtained.  It  should  be  added  that  for  nonlinear  cases 
there  is  no  guarantee  that  the  denominator  will  not  be  zero. 

The  importance  of  Q  in  Equation  (24)  is  to  specify  the  distance  of  the  current 
extrapolation  as  a  multiple  of  the  last  difference  vector,  xn+2  —  xn+i.  Though 
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Jennings’  two  suggestions  come  naturally  from  Aitken’s  formula,  values  for  Q  can 
be  found  by  other  effective  methods.  Chandler  (1987)  suggests  two  other  methods 
for  determining  Q,  that  based  on  results  obtained  on  test  problems  in  Chapter  XI 
work  as  well,  if  not  better  in  some  cases,  than  Jennings’  suggestions.  First,  define 
Q  as  the  quotient  of  Euclidean  norms: 


Q  =  (sign) 


!|Xn+2  ~  Xn+1  || 2 

||xn+2  —  2zn+i  +  Xnlh’ 


(26) 


where  sign  =  1  except  when  the  cosine  of  the  angle  between  the  difference  vectors, 
xn+2  ~  x„+i  and  xn+i  —  xn,  is  negative  or  the  norm  of  (x„+i  —  xn)  is  less  than  the 
norm  of  (xn+2  —  xn+i)-  In  these  cases  set  sign  =  —  1.  His  second  suggestion  is 

—  ||*n+2  —  Zn+1  ||  2 


Q  = 


(27) 


|xn+2  -  Xn+i||2  ±  ||Xn+l  -  *n||2’ 
where  the  denominator  is  a  sum  if  the  cosine  of  the  difference  vectors  is  negative 
and  a  subtraction,  otherwise. 

The  key  to  these  four  variations  is  that  they  all  extrapolate  along  the  vector 
%n+2  ~  Xn+ 1  •  Therefore,  it  is  essential  for  efficient  operation  of  the  algorithm  that 
the  cosines  of  the  difference  vectors  converge  to  either  plus  or  minus  unity.  If  the 
cosine  equals  plus  or  minus  unity,  then  all  four  methods  give  precisely  the  same 
results,  which  is  also  the  same  result  as  componentwise  Aitken.  However,  there 
exist  problems  where  all  four  suggested  methods  for  determining  Q  work  poorly. 
One  example  is  where  the  dominant  eigenvalues  axe  complex,  resulting  in  a  spiral¬ 
ing  convergence,  see  Figure  6.  A  second  example  is  where  a  single  real  eigenvalue 
does  not  dominate,  e.g.,  the  two-dimensional  problem  with  eigenvalues  of  0.8  and 
—0.8  where  the  cosine  of  the  angle  between  consecutive  difference  vectors  is  con¬ 
stant,  approximately  —0.42753,  see  Figure  10  (page  32).  However,  Jennings  pointed 
out  that  this  case  can  be  handled  very  well  by  taking  successive  pairs  of  the  ba¬ 
sic  iteration  to  form  new  basic  iterations  to  accelerate,  i.e.,  accelerate  the  sequence 
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Figure  10.  Trajectory  for  Eigenvalues  0.8,  —0.8 


{yn}  where  yn  =  x2n  or  yn  =  *2n+i>  n  =  0,1,...  This  effectively  squares  all  of  the 
eigenvalues  of  the  iteration  matrix  (Jennings,  1971).  Fortunately,  many  practical 
iterations  are  dominated  by  a  single  real  eigenvalue,  and  the  vector  Aitken  method 
with  the  Q  suggestions  of  Jennings  and  Chandler  are  often  effective  in  these  cases.  In 
addition,  Aitken’s  method  will  not  work  tremendously  well  on  divergent  sequences 
since  the  method  is  trying  to  approximate  the  solution  along  the  vector  xn+2  —  xn+1 . 
Chen  (1984)  suggests  interchanging  the  vectors  x„  and  xn+2  in  equation  (24).  Hence, 
the  extrapolation  will  be  along  the  vector  xn  —  xn+1 . 

Combining  Equation  (24)  with  Aitken’s  repeated  model,  a  semi-dynamic  al¬ 
gorithm  for  vectors  results.  Henceforth,  when  vector  Aitken’s  method  is  mentioned, 
it  is  referring  to  this  algorithm.  Unless  otherwise  stated,  Q  will  be  Jennings’  SDM. 
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AITKEN’S  SEMI-DYNAMIC  ALGORITHM  4.1: 

Find  a  solution  to  x  =  G(x)  given  the  initial  approximation  vector  xo.  Define 
the  terminology  “If  converged”  to  mean  “If  ||y  —  x||  <  Tol,  where  Tol  is  some 
predetermined  tolerance  value  and  x  and  y  are  the  vectors  of  the  current  step.” 

Step  1.  Compute  x\  =  G(x o)* 

Step  2.  Compute  x2  =  G(x i).  If  converged,  stop;  otherwise,  go  to  step  3. 

Step  3.  Find  y  =  x2  —  Q(x 2  —  xi),  where  Q  is  defined  by  Equation  (25),  (26), 
or  (27). 

Step  4.  Compute  z  =  G(y).  If  converged,  stop;  otherwise,  set  x0  =  y, 
x\  =  z  and  go  to  step  2. 

One  may  use  any  norm  desired  for  the  stopping  criterion.  There  are  advan¬ 
tages  and  disadvantages  for  all  of  them.  The  Euclidean  norm  will  show  a  smoothness 
in  the  differences  as  they  approach  the  tolerance  value.  However,  the  use  of  this 
norm  could  result  in  system  overflow  due  to  the  squaring  of  the  difference  compo¬ 
nents  unless  care  is  taken  in  computing  Q.  The  infinity  norm  is  simply  the  largest 
component  of  the  difference  vector,  y  —  x .  It  is  not  as  smooth  as  the  Euclidean 
norm  since  the  largest  component  may  be  a  different  component  from  one  iteration 
to  the  next.  However,  due  to  its  simplicity,  results  in  this  study  are  based  on  the 
infinity  norm  of  the  difference  vector. 

In  addition,  there  are  other  stopping  criteria  which  may  be  used: 

ll»-*.ll<(Tol)||y||,  j?*0,  and  ||F(y)||  <  Tol, 

where  F(x)  =  0.  Unfortunately,  difficulties  can  arise  no  matter  which  of  the  stop¬ 
ping  criteria  we  use.  For  example,  the  sequence  defined  by  xn  =  ^£=1(l/&)  is 
divergent  but  limn_00(xn  —  xn_i)  =  0.  This  sort  of  harmonic  convergence  is  “sub- 
linear”  (Brent,  1972)  and,  ordinarily,  is  never  encountered  in  practical  fixed  point 
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problems.  For  purpose  of  this  study,  the  stopping  criterion  for  all  algorithms  will 
be  as  stated  in  Algorithm  4.1. 

Numerical  Examples 


As  an  example  of  the  vector  Aitken  method,  consider  the  two-dimensional 
problem  of  finding  the  solution  to  the  system 

x  =  A(  sin  x)  +  B(  cosy)  and 

(28) 


y  =  A(cos  x)  —  2?(siny), 

where  A  =  0.7,  B  =  0.2,  and  x0  =  y0  =  0  (Henrici,  1964).  The  iteration  function 
then  becomes 


xn+1  =  A(sinxn)  +  B(cosyn)  and 
yn+ 1  =  A(cosxn)  —  B(sinyn). 

After  thirty-one  iterations,  we  obtain  the  solution  vector  to  five  decimal  places, 
(0.52652,0.50792).  If  we  apply  vector  Aitken  in  static  form,  we  obtain  the  solution 
to  the  same  degree  of  accuracy  in  12  iterations.  Table  5  (page  35)  shows  the 
Euclidean  norm  of  the  error  vector,  denoted  by  E ^  =  ||x^  —  s||2,  where  the  ith 
column  represents  the  sequence  derived  by  applying  Aitken’s  method  i  times.  The 
zero  column  is  the  original  iteration  sequence. 

Using  Aitken’s  method  semi-dynamically,  Algorithm  4.1,  the  same  results  are 
obtained  in  15  iterations,  Table  6  (page  36).  However,  Anderson  (1965,  p.  551)  says 
that  Aitken’s  method  “is  considerably  less  effective  if  applied  statically,  ...  than  it 
is  if  applied  dynamically.”  The  static  model  requires  32  extrapolations  to  only  6 
for  the  semi-dynamic  model.  Therefore,  the  benefit  of  3  fewer  iterations  is  loss  due 
to  the  time  required  to  perform  26  additional  extrapolations.  In  addition,  since 
the  number  of  column  sequences  are  not  known  beforehand  for  the  static  method, 
the  semi-dynamic  method  does  not  require  the  allocation  of  storage  space  to  ensure 
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TABLE  5 

AITKEN’S  STATIC  METHOD  APPLIED  TO  PROBLEM  (28) 
WITH  RESULTS  AS  EUCLIDEAN  NORM 
OF  ERROR  VECTOR 


n 

EW 

Ew 

Em 

£(M 

0 

0.53521 

1 

0.37883 

2 

0.23961 

0.25854 

3 

0.16527 

0.16557 

4 

0.11013 

0.07612 

0.01575 

5 

0.73253 

0.00578 

0.31493 

6 

0.04820 

0.00328 

0.00336 

0.15449 

7 

0.03156 

0.00154 

0.00182 

0.00181 

8 

0.02058 

0.00077 

0.00041 

0.00062 

0.00061 

9 

0.01338 

0.00032 

0.00005 

0.00015 

0.00062 

10 

0.00868 

0.00014 

0.00004 

0.00003 

0.00001 

0.00021 

11 

0.00563 

0.00006 

0.00001 

0.00001 

0.00001 

0.00001 

12 

0.00365 

0.00003 

0.00000 

0.00000 

enough  columns  for  convergence.  Thus,  the  semi-dynamic  model  is  usually  more 
efficient  than  the  static  model. 

Irons  and  Shrive  (1987)  made  a  modification  to  Aitken’s  method  for  scalars. 
Assume  that  we  have  the  two  relations  y\  =  G(y0)  and  y3  =  G(y3),  that  the  ratio 
of  consecutive  error  terms  is  constant,  and  that  s  is  the  limit  of  the  sequence  {t/n}. 
Then  the  following  is  true: 

_  Vi  ~  3  _  V3~  3 
VO  -  3  V2  ~  3 
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TABLE  6 

AITKEN’S  SEMI-DYNAMIC  METHOD 
APPLIED  TO  PROBLEM  (28) 


n 

En 

n 

En 

n 

En 

1 

0.37883 

6 

0.00393 

11 

0.00005 

2 

0.23961 

7 

0.00291 

12 

0.00004 

3 

0.17611 

8 

0.00189 

13 

0.00001 

4 

0.11835 

9 

0.00016 

14 

0.00001 

5 

0.00734 

10 

0.00010 

15 

0.00000 

Solving  for  s  gives 

_  _  (y»  -  -  jft) 

3  (yo  -yi)  -  (y2  -ya)' 

We  can  use  (29)  as  a  model  iteration  formula  for  estimating  the  limit  of  a 
sequence  for  which  the  ratios  of  consecutive  errors  converge  to  a  constant.  Given 
the  scalars  yo,yi,y2,  and  y3;  then 

(yn+2  —  yn+3)(yn+l  —  yn+3)  _  j 

y«+4  —  yn+3  -  7 - 7  7 - r  and 

(yn  yn+l)  (yn+2  yn+3) 

yn+s  =  <?(yn+4),  n  =  0,2,4,...  (30) 

Formula  (30)  gives  us  a  third  type  of  model,  a  fully  dynamic  method.  Study¬ 
ing  the  diagram  of  the  dynamic  model  in  Figure  11  (page  37),  we  see  that  each 
extrapolation  after  the  first  one  uses  only  one  additional  iterate  and  data  obtained 
from  previous  iterations.  A  dynamic  model  does  not  require  restarting  our  pro¬ 
cedure  by  generating  ail  necessary  iterates  from  the  latest  extrapolation.  We  will 
see  in  later  chapters  the  advantages  of  this  model.  Given  the  sequence  {x„},  we 
may  apply  Irons  and  Shrive’s  dynamic  model  after  two  iterations  by  setting  y0  =  x0, 
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yo  ->  yi 

yz  -*  V3 


EXT(jto,yi,y2,jte)  -*•  y4 


ys 

y4 


ys 

ys 


EXT(y2)y3,y4,y5) 


ye 


>  EXT(y4,ys,y6,y7)  ya  etc. 


y4  ys 

y«  -*■  y? 

Figure  11.  Diagram  of  Dynamic  Model 


yi  —  y2  =  Xi,  and  y3  =  x2.  Hence,  the  first  extrapolation  will  be  identical  to 
Aitken’s  first  extrapolation;  however,  the  equalities  stop  at  this  point. 

Consider  the  iterative  equation 

xn+1  =  2 sin(a;n),  n  =  0,l,...,  and  x0  =  1.  (31) 

Table  7  (page  38)  shows  a  comparison  of  the  results  obtained  when  Aitken’s  semi- 
dynamic  method  and  Irons  and  Shrive’s  fully  dynamic  method  are  applied  to  (31). 
The  column  headings  represent  the  current  iterated  value  (iter)  and  the  current 
extrapolated  value  (ext),  if  one  was  applied  on  that  iteration. 

Unsuccessful  Application 

In  Chapter  II,  it  was  noted  that  wrong  answers  can  sometimes  be  obtained 
when  trying  to  accelerate  a  sequence  by  Aitken’s  method.  Shanks  (1955)  found 
such  a  problem.  Let  {xn}  be  the  sequence  of  partial  sums  of  the  function 

f{l)  =  (l-*)(2-z) 

=  l  +  (3/2)z  +  (7/4)z2  +  (15/8)z3  +  ... 
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TABLE  7 


DYNAMIC  VS  SEMI-DYNAMIC  COMPARISON 
FOR  PROBLEM  (31) 


Semi-Dyi 

namic 

Dynamic 

n 

iter 

ext 

iter 

ext 

0 

1.000000 

1.000000 

1 

1.682942 

1.682942 

2 

1.987436. 

1.915372 

1.987436 

1.915372 

3 

1.882438 

1.882438 

1.892686 

4 

1.903662 

1.895344 

1.897278 

1.895462 

5 

1.895590 

1.895514 

1.895494 

6 

1.895434 

1.895495 

7 

1.895494 

If  z  =  4,  the  series  diverges.  However,  it  was  shown  in  Chapter  II  that  Aitken’s 
method  can  obtain  the  antilimit  of  a  divergent  sequence.  Table  8  (page  39)  shows 
the  static  results  when  Aitken’s  method  is  applied  to  this  problem.  Column  i  is 
the  sequence  obtained  by  the  ith  application  of  the  method.  We  see  that  the  later 
columns  converge  to  7/27  =  0.25926...  However,  the  value  of  /(4)  is  1/3.  Thus, 
Aitken’s  method  did  converge,  but  to  the  “wrong”  value.  Shanks  showed  that 
erroneous  results  in  this  problem  will  be  obtained  only  for  z  =  4.  Hence,  one  may 
feel  confident  that  wrong  results  are  few  and  far  between  in  practical  applications, 
but  do  exist.  However,  Lubkin  (1952)  did  prove  that  if  any  two  consecutive  columns 
of  the  Aitken’s  table  both  converge,  then  they  converge  to  the  same  limit.  Hence, 
if  all  columns  converge,  then  they  all  converge  to  the  correct  limit.. 
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TABLE  8 

CONVERGENCE  TO  WRONG  ANSWER  BY 
AITKEN’S  STATIC  METHOD  APPLIED 
TO  SHANKS’  EXAMPLE 


n 

0 

1 

2 

3 

4 

0 

0 

1 

1 

2 

7 

-0.2000 

3 

35 

-0.6364 

4 

155 

-1.5217 

0.2241 

5 

651 

-3.2979 

0.2437 

6 

2667 

-6.8526 

0.2519 

0.2589 

7 

10,795 

-13.9634 

0.2557 

0.2589 

8 

43,435 

-28.1854 

0.2575 

0.2592 

0.2593 

CHAPTER  V 


SHANKS’  TRANSFORMATIONS  FOR  SCALARS 

In  this  chapter  the  general  framework  for  deriving  the  remaining  acceleration 
techniques  will  be  established.  The  motivation  will  be  similar  to  Shanks’  (1955) 
development  of  his  e*.  transformation  for  scalar  sequences.  First,  consider  a  variety 
of  typical  scalar  sequences  {zn}:  convergent,  divergent,  monotonic  and  oscillatory. 
Plotting  the  sequence  elements  versus  n  and  connecting  them  with  a  smooth  curve 
give  graphs  similar  to  the  samples  shown  in  Figure  12. 


♦si 


♦ 


Figure  12.  Plots  of  Typical  Sequences 
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Shanks  observed  that  the  graphs  all  look  like  the  graphs  of  what  he  called 
“physical  transients.”  By  this  term  he  meant  a  physical  quantity,  p,  as  a  function 
of  time  in  the  form: 

k 

p(t)  =  s  +  ebit, 

i= 1 

where  the  Vs  axe  complex  numbers,  bi  ^  0.  He  referred  to  them  as  “transients” 
so  that  he  could  apply  the  term  in  a  more  general  sense  for  both  convergent  and 
divergent  sequences.  Since  p(t)  is  an  exponential  function  with  6;  an  element  of 
the  complex  numbers,  the  set  of  functions  also  include  the  trigonometric  functions. 
Shanks  represented  the  sequences  {z„}  as  if  they  were  “mathematical  transients,” 
a  function  of  n  in  the  form 


k 

xn  =  s +  ^04  q?,  qi^  0,1, 

i=l 

where  s,<Zi,  and  qi  are  constants  independent  of  n  and  qi  ^  qj  for  i  /  j.  There¬ 
fore,  his  “mathematical  transient”  equation  is  identical  to  the  relationship  (13)  for 
scalars.  Once  again,  we  have  the  concept  that  s  is  either  the  limit  or  the  antilimit 
of  {«„}  depending  upon  the  moduli  of  the  V s.  As  stated  earlier,  Shanks  referred 
to  a  divergent  sequence  as  “diverging  from  s”  (1955,  p.  7). 

Shanks’  proposed  method  of  approximating  s  was  to  solve  the  2k  +  1  system 
of  nonlinear  equations, 


k 

xr  =  Bh<n  +  Y,  (H  qri ,  n  <  r  <  n  +  2k, 

i=  1 

for  Bk,n  with  a,-,  qi,  i  =  1, . . . ,  k,  the  rest  of  the  unknown  values.  Here  „  is  taken 
to  be  an  approximation  of  s.  Shanks  determined  that  the  solution  Bk,n  could  be 
represented  in  the  determinant  form 
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Xn 

zn+l 

*  *  * 

Xn+k 

Un 

U„+l 

.  .  . 

Un+k 

Un+k-1 

Un+k 

Un+2k-l 

1 

1 

1 

Un 

Un- (-1 

.  .  . 

Un+k 

Un+k-\ 

Un+k 

Un+2k-l 

(32) 


n  =  0, 1, . . and  where  un  =  xn+i  —  xn.  Therefore,  Shanks  derived  a  new  sequence 
{B*  n}  defined  by  (32)  where  k  is  a  nonnegative  integer  and  for  which  the  denom¬ 
inator  does  not  vanish.  If  the  denominator  vanishes  for  n  =  m  and  the  numerator 
does  not,  then  Bk>m  is  assigned  the  value  oo.  If  both  numerator  and  denominator 
vanished  for  n  =  m,  then  Bk<m  —  Hfc-i.m •  He  wrote  the  transforms  in  operator  form 
as 


ek(xn)  =  Bk,n 

Shanks  called  ek  “the  fc’th  order  transform  of  {An},”  (Shanks,  1955,  p.  2) 
where  {.An}  =  {xn}.  There  are  two  transforms  that  need  to  be  identified.  The  first 
one  is  k  =  0  where  e0(xn)  =  xn-  For  the  second  one,  letting  fc  =  1  we  have 


Xn 

*n+l 

Un 

Un+ 1 

®n(®n+2  ®n+l)  ®n+l  ( "Eri-l-l  Xn) 

1 

1 

(a-n+2  S-n+l  )  (-^n+l  Xn ) 

Un 

Un+1 

XnXn+2  Xn+I 
Xn+ 2  2ln+i  -f-  Xn 


n  =  0, 1, . . . 


(33) 
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Comparing  (33)  with  (22)  we  see  that  the  two  right-hand  expressions  are  equivalent. 
Therefore,  the  first  extrapolated  sequence  of  Aitken’s  A2  method  and  Shanks’  ex 
transformation  are  equivalent  for  scalars.  However,  this  is  not  true  for  a  second, 
third,  or  kth  application  of  the  methods.  In  the  next  chapter,  it  will  be  shown  how 
the  sequences  of  scalars  obtained  by  higher  order  applications  (second,  third,  etc.) 
of  Aitken’s  method  can  be  obtained  from  the  e&  transformations.  However,  the 
e/,  transformation  sequences  will  be  derived  by  a  technique  different  from  Shanks’ 
determinant  method.  As  one  can  clearly  see,  the  larger  the  value  of  k ,  the  more 
complicated  the  solving  for  Bk,n,  since  two  (k  -f  1)  x  (k  +  1)  determinants  must 
be  computed.  Hence,  a  method  of  obtaining  similar  results  without  the  use  of 
determinants  would  be  a  valuable  asset.  In  the  next  chapter,  such  a  technique  is 
presented. 


CHAPTER  VI 


WYNN’S  EPSILON  AND  MODIFIED  EPSILON  METHODS 
Theory  for  Scalar  Epsilon  Method 

Wynn  (1956)  derived  the  epsilon  (e)  algorithm  to  accelerate  a  sequence  of 
scalars.  His  simple  algorithm  effected  Shanks’  ek(xn)  transformation  without  the 
use  of  determinants.  Wynn  showed  that  he  could  calculate  ei(®„)  directly  from  the 
elements  of  {xn}  and  e^(xn)  directly  from  the  ei_i(xn)  elements  and  values  deter¬ 
mined  from  the  i(xn)  elements,  i  ~  2,3,...  He  used  the  symbol  to  represent 
his  new  values,  where  the  k  subscript  refers  to  a  column  number  and  the  n  super¬ 
script  refers  to  a  diagonal  number.  His  method  will  be  derived  for  both  scalars  and 
vectors. 

For  the  scalar  case,  the  values  e £  are  determined  from  a  given  sequence  by 
setting  the  initial  two  column  values  as 

=  0,  £q  =  xn,  n  =  0,1,..., 

and  by  the  relationship 


The  quantities  e*  may  be  arranged  as  shown  in  Figure  13  (page  45).  Note  that  the 
four  quantities  of  (34)  are  located  at  the  four  corners  of  a  lozenge,  as  indicated  for 
n  =  0,k  =  2  in  Figure  13.  Therefore,  a  quick  way  of  remembering  how  to  find  the 
right  side  entry  of  the  lozenge  is 

Right  =  Left  +  (Bottom  —  Top)-1  . 
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£°o  =  xQ 
ej  =  xi 
e\  =  x2 

So  =  *3 

Sq  =  z4 


£° 

£2 


e\ 


£ 


0 

4 


Figure  13.  Wynn’s  Epsilon  Arrangement 


The  odd  and  even  numbered  columns  Me  quite  varied  in  the  information 
they  give.  The  odd  subscript  columns  normally  diverge  and  give  no  directly  useful 
information  as  to  the  limit  or  antilimit  of  the  sequence.  However,  one  can  obviously 
see  that  they  are  vital  as  they  are  used  to  determine  the  next  even  column.  The 
even  subscript  columns  will  often  converge  to  the  desired  limit  or  antilimit  of  {xn} 
and  will  do  so  more  quickly  than  the  original  sequence.  However,  the  key  sequence 
for  convergence  is  the  diagonal  sequence  whose  elements  aure  4n>  to  =  0, 1, . . .  This 
sequence  will  most  often  converge  not  only  the  quickest,  but  in  some  cases,  will 
converge  even  when  each  even  numbered  column  sequence  diverges. 

In  order  to  find  eJ|,  one  must  first  have  computed  £5,  £5,  and  e\.  This  even¬ 
tually  leads  to  the  fact  that  before  £g  can  be  found,  all  elements  in  the  first  six 
“cross-diagonals”  must  be  determined.  The  term  cross-diagonals  refers  to  the  di¬ 
agonals  of  Wynn’s  epsilon  arrangement,  Figure  13,  which  rise  as  one  moves  from 
left  to  right  in  the  figure  and  where  n  +  k  is  constant.  Wynn  (1964)  suggested 
that  to  prevent  the  use  of  unnecessary  storage,  computation  of  the  elements  can  be 


46 


made  by  computing  one  cross-diagonal  at  a  time,  using  only  the  data  saved  from 
the  previous  cross-diagonal. 

Therefore,  to  use  his  technique  as  an  acceleration  method,  elements  of  the 
arrangement  diagram  are  computed  one  cross-diagonal  at  a  time  until  the  e02m  ele¬ 
ment  has  the  desired  precision  of  accuracy  as  measured  by  ||G(£°m)  —  £2mll  <  Tol, 
where  G  is  the  iteration  function.  Wynn’s  theorem  relating  £fc(x„)  (or  e£)  to  efe(xn) 
follows. 

THEOREM  6.1:  If  £2m(xn)  =  em(z„)  and  £2m+i(zm)  =  (em(un))"\  where 

Un  —  Xn+l  then  £jfl(jn)  =  £*— 1  (*n+l  )  4"  (£j(Xn+l  )  £*(xn))  )  •S  —  1, 2, . . . 

Wynn’s  proof  was  by  mathematical  induction  (Wynn,  1956,  p.  92-94).  He 
proved  the  equality  by  showing  that 

£«+i(xn)  -  £.-i(xn+i)  and  (s,(xn+1)  -  e,(xn))_1 

axe  equivalent  expressions  in  determinant  form.  Hence,  in  certain  cases,  the  se¬ 
quence  £°m,  771  =  0,1,...,  converges  to  the  limit  or  antilimit  of  {xn}  and  the  con¬ 
vergence  is  more  rapid  than  the  original  sequence. 

For  an  example  of  how  Wynn’s  epsilon  method  works,  consider  Equation  (7). 
Table  9  (page  47)  shows  the  epsilon  arrangement  for  five  iterations.  Hence,  Wynn’s 
£  method  converges  in  the  same  number  of  iterations  as  Aitken’s  static  method.  In 
addition,  this  example  will  be  used  to  illustrate  how  the  scalar  sequences  obtained 
by  more  than  one  application  of  Aitken’s  A2  method  can  be  computed  by  the  use  of 
Wynn’s  epsilon  method.  Once  column  two  has  been  computed,  calculate  the  next 
two  columns  as  if  the  second  column  were  column  zero.  In  other  words,  when  com¬ 
puting  column  three,  treat  column  one  as  if  every  element  were  zero.  Repeat  this 
process  with  columns  four,  five,  and  six,  etc.  When  all  is  done,  the  even  numbered 
columns  will  be  the  same  sequences  as  would  be  derived  by  repeated  applications 
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TABLE  9 


WYNN’S  EPSILON  ARRANGEMENT  FOR  EQUATION  (7) 


71 

eS 

£1 

£2 

c3 

Fn 

£4 

0 

0.5000000 

3.586790 

1 

0.7788008 

-9.867936 

0.7044777 

-1141.7298 

2 

0.6774630 

28.400377 

0.7035942 

-8964.4053 

0.7034663 

3 

0.7126738 

-80.404595 

0.7034830 

-73073.1050 

0.7034674 

4 

0.7002567 

228.937720 

0.7034693 

5 

0.7046047 

of  Aitken’s  A2  method.  Table  10  (page  48)  shows  the  results  of  this  procedure  for 
columns  two,  three,  and  four.  One  may  check  that  the  columns  labeled  ejj'  and  e£' 
are  identical  to  columns  two  and  three  of  Table  3,  which  were  obtained  by  applying 
Aitken’s  method  to  this  same  problem. 

Theory  for  Vector  Epsilon  Method 

Now  consider  {xn}  as  a  sequence  of  vectors.  Before  Wynn’s  algorithm  can  be 
applied  to  a  vector  sequence,  the  inverse  of  the  vector  x  =  (xi,x2, . . .  ,  xm)r  must 
be  defined.  Wynn  (1962)  discusses  two  possible  inverses: 

(1)  Primitive  Inverse: 

r 

x~x  =  (x^jX,1,...,®-1)  ,  where  x,- ^  0  for  all  i. 

For  Xi  =  0,  take  x,-1  =  0. 

(2)  The  Samelson  Inverse: 

m 

X  ~  (*ll  •  l  *m)  » 

r=l 
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TABLE  10 

AITKEN’S  A2  SEQUENCES  DERIVED 
FROM  WYNN’S  ARRANGEMENT 


n 

Pn> 

£n' 

£1 

£n‘ 

£2 

1 

0.7044777 

-1131.8619 

2 

0.7035942 

-8992.8058 

0.7Q34669 

3 

0.7034830 

-72992.7007 

0.7034674 

4 

0.7034693 

where  xr  is  the  complex  conjugate  of  xr  and  x  is  not  the  zero  vector.  Define  the 
zero  vector  as  the  inverse  of  itself.  Samelson’s  inverse  is  equivalent  to  the  Moore- 
Penrose  generalized  inverse  of  x  considered  as  a  m  x  1  matrix.  Hence,  it  will  always 
be  referred  to  as  such.  The  primitive  inverse  ignores  the  relationship  between  the 
scalar  sequences  of  different  components.  In  addition,  it  will  frequently  have  major 
problems  in  that  one  or  more  of  the  reciprocals  will  be  quite  large  numerically  due 
to  a  denominator  very  close  to  zero.  Further  more,  as  far  as  is  known,  the  primitive 
inverse  seldom  gives  better  results  than  the  generalized  inverse  (Smith,  Ford,  Sidi; 
1987,  p.  223).  Therefore,  the  generalized  inverse  is  more  useful  and,  hence,  all  work 
in  this  study  involving  inverses  of  vectors  will  use  generalized  inverses. 

WYNN’S  VECTOR  EPSILON  ALGORITHM  6.2: 

Given  the  iteration  equation  xn+1  =  G(xn)  and  the  initial  vector  x0.  Define 
i"h  =0,  n  =  1, 2, . . . ,  f  o  =  xo,  and  set  n  =  1. 

Step  1.  e  J  =  xn. 

Step  2.  For  k  =  0  to  n  -  1,  find  =  e*£-i  +  -  £ 1  • 
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Step  3.  If  n  is  even  and  ||e*°  —  £**_2||  <  Tol  or  if  n  is  odd  and 

—  £**_3||  <  Tol,  then  go  to  step  4;  otherwise,  n  =  n  4- 1  and 
go  to  step  1. 

Step  4.  Find  y  =  G(£*°)  if  n  is  even  or  y  =  <j(£*^_1)  if  n  is  odd.  If 
||y  —  e  °  (or  £**_x)||  <  Tol,  then  stop;  otherwise,  n=n+l  and 
go  to  step  1. 

Through  the  work  of  Cheng  and  Hafez  (1959),  the  epsilon  method  can  be 
modified  to  make  a  semi-dynamic  model.  Using  only  the  initial  vector  and  the  first 
two  iterates,  the  first  extrapolated  term,  e  is  found  by  Equation  (34).  Using 
this  vector  as  a  new  initial  vector,  two  new  iterates  are  generated  and  another 
extrapolated  vector  determined.  This  pattern  is  continued  until  an  iteration  has 
the  desired  precision  of  accuracy  as  measured  by  ||G(e£)  —  £  £||  <  Tol.  A  diagram 
of  the  model  is  shown  in  Figure  14  (page  50).  This  semi-dynamic  model  of  Wynn’s 
epsilon  method  is  called  the  Modified  Epsilon  Method. 

MODIFIED  VECTOR  EPSILON  ALGORITHM  6.3: 

Given  the  iteration  equation  xn-i  =  G(xn)  and  the  initial  vector  xo.  Define 

e "i  =  0,  n  =  0,1,2;  £*□  =  x0;  and  el  =  G(x0)- 
Step  1.  Compute  e\  =  G(e  q). 

Step  2.  Compute  e  ?,  e  J,  and  s  2  by  equation  (34). 

Step  3.  Compute  =  ^(^*2)’  ^  11^*2  —  ^*2!!  <  Tol,  then  stop;  otherwise, 
set  =  £*2,  el  =  £*2,  and  go  to  step  1. 

The  modified  epsilon  method  as  described  in  Algorithm  6.3  and  shown  in 
Figure  14  is  just  one  version  of  several.  The  method  shown  is  of  order  one  with 
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e\ 

e  0  —  e 

e  2  ~  e  o 

£  ,c 
e  i 

e\ 

e  o 

£  fO  _  -  //0 
6  2  ~  e  0 

s  i 

etc. 

Figure  14.  Diagram  for  Modified  Epsilon  Method 


respect  to  the  initial  value  =  x0.  Order  two  with  respect  to  uses  the  vectors 
through  e  o  to  determine  and  then  sets  =  e°  until  convergence.  Higher  order 
methods  continue  in  a  similar  fashion;  however,  orders  less  than  four  seem  to  work 
the  best,  especially  for  divergent  sequence  to  keep  the  iteration  from  diverging  too 
quickly.  As  was  the  case  with  Aitken’s  method,  the  semi-dynamic  model  is  more 
efficient  than  the  original  cross- diagonal  static  model  due  to  the  time  and  storage 
required  to  compute  the  triangular  arrangement  for  problems  of  large  dimension. 

Other  variations  start  with  a  different  element  of  the  original  sequence  as 
the  initial  value.  Here,  a  desired  number  of  elements,  say  i,  of  the  sequence  are 
skipped  and  the  process  begins  with  fj,.  This  procedure  is  useful  if  the  multiplicity 
of  the  eigenvalue  zero  is  known,  even  though  exact  zero  eigenvalues  are  rare  in 
practical  application.  Skipping  the  number  of  elements  equal  to  the  multiplicity 
of  the  eigenvalue  zero  will  result  in  faster  convergence.  Another  more  practical 
purpose  for  skipping  a  set  number  of  iterates  is  when  under-relaxation  is  used  to 
move  various  eigenvalues  closer  to  zero. 
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Numerical  Examples 

Consider  this  original  problem  of  Equation  (4)  where 

1.0  0.1  „  _ 

A=  and  b  =  (1.2,  -2.0)r.  (35) 

-0.5  0.4 

-  m 

The  solution  to  this  problem  is  (10.4,  —12.0).  Using  the  iterative  Equation  (5) 
and  the  initial  vector  (1,1),  Wynn’s  vector  epsilon  method  determines  the  solution 
with  only  four  iterations.  Table  11  (page  52)  shows  the  Euclidean  norm  of  the 
error  vectors  for  the  even  numbered  columns.  Remember  that  the  even  numbered 
columns  in  Wynn’s  arrangement  are  the  only  valid  sequences.  By  checking  the  rate 
of  convergence  of  the  first  column,  one  can  easily  see  that  Wynn’s  method  took  a 
very  slow  converging  sequence  and  accelerated  it  tremendously.  Table  12  (page  52) 
shows  the  results  of  applying  the  semi-dynamic  (modified)  model  to  problem  (35). 
Results  are  shown  for  orders  one,  two,  and  three. 

Wynn  conjectured  and  McLeod  (1971),  Theorem  6.4,  and  Gekeler  (1972), 
Theorem  6.5,  proved  the  following  theorems  which  were  thoroughly  discussed  by 
Brezinski  (1974). 

THEOREM  6.4:  Let  the  relation 

k  h 

53  ^  xn+i  =  i*53a»>  n  =  0,1,... 

«'=o  »=o 

hold  for  the  initial  values,  where  the  coefficients  a*  axe  real,  a*  ^  0,  s  and  xn  axe  m- 
dimensional  vectors  over  the  complex  numbers,  and  the  vectors  e*"  are  determined 
by  (34)  and  exist  for  n  +  r  <  2k.  Then 

k 

£*2 k  =  *  for  every  n  if  53  a*  i1  0,  and 

«=o 

£2 k  =  0  for  every  n  otherwise. 


TABLE  11 


EUCLIDEAN  NORMS  OF  WYNN’S  EVEN  NUMBERED 
COLUMN  ERROR  VECTORS  FOR  PROBLEM  (35) 


n 

Eno 

E* 

0 

16.042443 

1 

12.791403 

7.066314 

2 

10.710378 

5.960512 

0.000000 

3 

9.245948 

4.807537 

4 

8.127063 

TABLE  12 

EUCLIDEAN  NORMS  OF  ERROR  VECTORS  FOR 
MODIFIED  £  METHOD  FOR 

PROBLEM  (35) 

n 

Euclidean  Norms  for  Order 

1 

2 

3 

0 

16.042443 

16.042443 

16.042443 

1 

12.791403 

12.791403 

12.791403 

2 

10.710378 

10.710378 

10.710378 

3 

6.127186 

9.245948 

9.245948 

4 

5.399512 

8.127063 

8.127063 

5 

2.481779 

0.000000 

7.217934 

6 

2.078020 

6.448139 

7 

1.188792 

0.000000 

53 


THEOREM  6.5:  If  the  vector  ^-algorithm  is  applied  to  vectors  produced  by  the 
linear  system  (4)  where  A  is  a  real  matrix  such  that  (I  —  A)  is  nonsingular,  then 

e°2k  =  (I  -  A)-'b  =  s,  (36) 


where  k  is  the  degree  of  the  monic  minimal  polynomial  of  A  with  respect  to  the 
vector  Xq  —  s ;  that  is,  k  is  the  smallest  integer  such  that  there  exists  the  polynomial 

k 

p{y )  =  y\  Pk  =  i 

»=0 


where 


p(A.)(i0  -  s)  =  0. 

Equation  (36)  can  be  generalized  to 


,=*n+4 


where  0  <  q  <  r  for  r  equal  to  the  multiplicity  of  the  eigenvalue  zero  of  the  matrix 
A  (Brezinsld,  1974). 

The  significance  of  this  fact  will  primarily  be  seen  in  later  chapters  as  we  look 
at  other  acceleration  techniques.  However,  the  results  of  Problem  (35)  shown  in 
Table  11  illustrate  this  principle.  It  cam  be  shown  that 


p{y)  =  y2  ~  l-4y  +  0.45 

is  the  minimal  polynomial  of  A  with  respect  to  Xo  —  s.  Hence,  k  =  2  and,  there¬ 
fore,  e  °  should  be  equal  to  s  if  rounding  errors  are  not  considered.  The  minimal 
polynomial  of  the  matrix  A  will  also  be  considered  more  in  later  chapters. 


CHAPTER  VII 


THE  MINIMAL  POLYNOMIAL  EXTRAPOLATION  METHOD 

Theoretical  Aspect 

From  this  point  on  the  focus  of  this  study  will  center  on  vector  sequences  only. 
The  first  method  to  consider  is  a  method  developed  by  Cabay  and  Jackson  (1976). 
They  derived  a  polynomial  extrapolation  method  for  finding  the  limit  (antilimit) 
of  a  vector  sequence  {af}  governed  by  the  linear  iteration  (5).  They  assume  that 
(7  —  A)-1  does  exist  so  that  the  limit  is  the  unique  solution  of  equation  (4). 

The  key  item  in  their  method  is  the  minimal  polynomial  p(y)  of  A  which 
annihilates  u0  —  X\  —  af0;  in  other  words,  p  is  that  unique  monk  polynomial  of  least 
degree  such  that 

p(A)u o  =  0.  (37) 

Therefore,  their  technique  is  referred  to  as  the  minimal  polynomial  extrapolation 
(MPE).  To  derive  their  algorithm,  let  i*be  the  solution  of  (4)  and  define 

u  =  3  —  x0  and  un  =  ®n+i  —  xn.  (38) 

Hence, 

un+i  =  Au„,  (7  -  A)u  =  u0,  and 

(39) 

(I-A)(f-  xj)  =  Uj. 

Let 

sp(v)  =  E|  E  ci )  y'»  where 

*=0  \J=*+ 1  / 

k 

p{y )  =  Y,cjy3'  c*>  =  i,  (40) 

j=0 
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Hence, 


(/  -  A)Sp(A) 


[(ci  +  .  .  .  +  C]t)I  +  (c2  +  .  .  .  +  Ck)A  +  .  .  .  +  CkAk  l]  — 

{(ci  4- ...  4-  Ck)AA°  4-  (02  +  . . .  +  Ck^AA1  4-  . . .  +  CkAA  1 } 
(ci  4-  C2  4"  •  •  •  4-  £k)I  —  (ci  4-  C2  4-  •  •  •  4-  Ck)Ak 


=  P{I)~P(A). 


Assuming  p(A)  annihilates  u0,  then 


(J  -  A)SP(A) uo  =  ( p(I )  -  p(A))u0  =  p(I)u0. 

Using  equations  (38)  and  (39)  and  simplifying  give 

U  =  3-  xQ  =  (I  -  A)-1u0  =  p(l)-15p(A)zf0- 
However,  it  is  also  true  that 

Ui  =  xi+1  —  Xi  =  A£i  —  Axi^i  =  ...  =  A'u0. 
Hence,  (41)  becomes 

3-x0  =  p(l)'1 


fc-i  (  k  \ 

EE  c)  j  A‘ 

,«= o  \j=i+ 1  / 


u0 


=  P(l) 


-l 


fc-i  /  k 

E  E  <*|* 

L»=o  \;=»+l 


Therefore,  solving  for  T  gives 


(41) 


(42) 


3  =  X0  + 


k- 1  . 

(  k  \ 

E| 

E  <* 

|  Ui 

_»=° 

\j=i+i  ) 

m 

/?(!) 


(43) 


Since  (/  —  A)-1  exists,  A  has  no  eigenvalue  at  unity  and  thus  p(l)  ^  0.  There¬ 
fore  u  can  be  found  after  &  4-  1  iterations,  provided  the  annihilating  pol  nial  ex¬ 
ists.  Cabay  and  Jackson  made  no  attempt  to  produce  the  minimal  polynonual  p(y). 
Instead,  they  found  an  almost-annihilating  polynomial  a(y)  =  °>  y\  <*(1)  ^  0 


and  ay  =  1,  such  that  YJlLo  g»  for  6  relatively  small.  Once  the  a,’ s  are  found, 

the  extrapolated  vector  is  determined  by  calculating  u  and  adding  the  result  to  x0- 
One  method  for  solving  the  a;’ s  is  to  minimize  the  norm  by  using  a  least  squares 
technique.  Hence,  we  solve 

fc'-i 

^2  CLi  Ui  =  —Uy. 

*=0 

Another  approach  (Sidi,  1986,  and  Sidi,  Ford,  and  Smith,  1986)  in  developing 
the  MPE  method  is  to  begin  with  the  minimal  polynomial,  p(y ),  of  A  with  respect 
to  u0,  Equation  (40).  Using  (37)  and  (42), 

k  k 

0  =  p(A)u0  =  ^2  Cj  A:  u0  =  Yl  ci  “j- 

j—0  j= 0 

So  the  unknown  coefficients  of  the  polynomial  p  are  =  1  and  the  components  of 
the  vector  c  =  (c0,  Ci, . . . ,  Cfc_i)r  which  solves  the  system  of  equations 

Uc  =  — tffe,  (44) 

where  U  is  the  m  x  k  matrix  defined  by 

U  =  [u0,ui, . . .  ,Ufc-i]. 

If  k  <  m,  then  there  are  more  equations  in  the  system  than  unknowns;  however,  we 
have  shown  consistency.  Therefore,  the  unique  solution  can  be  found. 

We  now  express  any  element  Xj  in  terms  of  x0  by  using  the  fact  that  if  i*  is 
the  solution  of  (4)  then  s  =  (I  —  A)-1  b.  So 

Xj  =  A3  x o  +  (/  +  A  -+-...  -+•  A3  1 )b 

=  A3x o  +  (f  -  A3){I  -  A)~lb  =  A3x o  +  (/  -  -4J)i* 

=  A3(x o  -  s')  +  i*. 

Hence,  A3(x0  —  s  )  =  xj  —  s. 
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Smith,  Ford,  and  Sidi  (1987)  showed  that  the  minimal  polynomial  of  A  with 
respect  to  the  vector  xj  —  sis  the  same  polynomial  as  that  for  Uj,  for  every  j,  and 
thus  true  for  j  =  0.  (Hence,  p(y)  is  the  same  minimal  polynomial  that  was  discussed 
in  Chapter  VI  when  it  was  shown  that  e^k  =  s  for  k  the  degree  of  p(y).)  So 

0  =  Yjci  Aj  “0  =  1] °j  AJ  (xQ  -  s)  =  J] Cj  (xj  -  s  ) 

j= 0  j= 0  j=0 

=  Ec>*-mEc;|  •  (46) 

j= 0  \j=o  ) 

Since  unity  is  assumed  not  to  be  an  eigenvalue  of  A,  =  p(l)  ^  0-  Hence,  s 

is  computed  directly  from  (46). 

The  above  proof  can  also  be  shown  for  the  starting  vector  xn  instead  of  x0- 
Therefore,  a  theorem  for  any  k  +  1  consecutive  terms  of  a  sequence  was  proven  by 
Smith,  Ford,  and  Sidi  (1987). 

THEOREM  7.1:  For  any  k+ 1  consecutive  terms  of  the  sequence  {x},  say  xn,  xn+1 , . . . 
xn+fe,  we  have 

Ec>  *"+i  =  *  (Ecl)  •  (47) 

:= o  \i=o  / 

where  Cj,j  =  0, . . . ,  k,  are  defined  by  equation  (40). 

If  (43)  is  rewritten  in  terms  of  the  x/s,  we  have 

*  =  E  li  where  h  =  ci!  fE  c’V 

j= o  \»=o  / 

If  r  is  the  multiplicity  of  the  eigenvalue  zero,  then  r  terms  on  each  side  of  (47) 
are  zero.  Therefore,  if  r  is  known  or  suspected  to  be  positive,  there  is  an  advantage 
of  starting  the  k  4-1  consecutive  terms  at  xn,n  >  0,  instead  of  Xo-  Preferably,  we 
should  start  at  xr. 

To  this  point,  the  discussion  has  dealt  with  finding  the  solution  of  a  linear  sys¬ 
tem.  For  linear  problems  of  large  dimension,  the  degree  of  the  minimal  polynomial 
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may  be  difficult  to  determine.  For  nonlinear  problems,  the  annihilating  polynomial 
changes  for  each  iteration  and  the  limit  cannot  be  obtained  in  a  finite  number  of 
iterations.  Therefore,  a  small  value  of  k  is  chosen  and  Equation  (43)  is  used  as 
a  model  for  approximating  the  solution.  Repeating  the  process  until  convergence 
with  each  new  initial  vector  set  equal  to  the  last  computed  extrapolated  vector,  a 
semi-dynamic  method  is  developed  for  solving  large  linear  problems  and  nonlinear 
problems.  Even  for  small  linear  problems  where  the  exact  k  is  known,  rounding 
errors  may  prevent  obtaining  the  solution  to  the  desired  accuracy  on  the  first  ex¬ 
trapolation.  Hence,  the  semi-dynamic  model  is  used  for  all  types  of  problems. 

MINIMAL  POLYNOMIAL  EXTRAPOLATION  (MPE)  ALGORITHM  7.2: 

Given  the  sequence  xn+i  =  G(xn),  the  initial  value  x0,  and  the  positive  integer  k. 

Step  1.  Generate  a?i,  x2,  •  •  • ,  *fc+i  by  the  function  G. 

Step  2.  Compute  U  and  tf*  by  use  of  (45)  and  (38). 

Step  3.  Compute  c  from  (44)  and  set  c&  =  1. 

Step  4.  Compute  s  from  (47)  where  xn+;  =  Xj. 

Step  5.  Generate  y  =  G(s).  If  ||y  —  i*||  <  Tol,  then  stop;  otherwise,  set 
x0  =  3  and  go  to  step  1. 

Theoretical  Application  to  Numerical  Problems 

Let  us  consider  Problem  (35)  again.  In  Chapter  VI,  it  was  shown  that  Wynn’s 
vector  epsilon  method  obtains  the  solution  in  four  (2»c)  iterations  since  the  degree 
of  the  monic  minimal  polynomial  is  k  =  2.  Since  this  example  is  a  linear  problem, 
according  to  theory,  the  MPE  method  should  compute  the  solution  in  k  +  1  =3 
iterations.  This  is  the  case  as  the  solution,  (10.4,-12.0),  is  found  to  six  decimal 
place  accuracy  after  only  three  iterations  and  one  extrapolation.  Therefore,  in 
theory,  whereas  Wynn’s  Vector  Epsilon  method  requires  2k  iterations  to  obtain  $, 
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the  MPE  method  can  obtain  it  in  only  k  +  1  iterations.  If  k  is  large,  one  can  see  a 
major  advantage  of  the  MPE  method.  However,  the  larger  the  value  of  k ,  the  larger 
the  dimension  of  our  matrix  U  in  determining  the  coefficient  vector,  c.  If  the  MPE 
method  is  applied  to  (35)  with  k  =  1,  the  extrapolation  procedure  still  converges; 
however,  it  takes  78  iterations  to  determine  the  solution  to  six  decimal  places. 

For  larger  dimensional  problems,  k  usually  should  be  chosen  such  that 
2  <  k  <  5.  There  are  two  reasons  for  this  restriction  of  k.  The  first  reason  is  the 
storage  space  and  the  time  required  in  working  with  large  dimensional  matrices.  The 
second  reason  is  given  by  Anderson  (1965,  p.  555),  “the  power  of  an  iterative  method 
increases  slowly  with  degree  for  M  >  3  since  the  “early,”  poor  approximations  are 
not  samples  of  significant  information  content  ...”  (Here  Anderson’s  M  refers  to  k). 
Results  found  from  the  numerical  test  problems  in  Chapter  XI  show  that  no  one  k 
is  the  best  choice  for  all  problems. 

Variations  for  Convergent/Divergent  Sequences 

Since  Equation  (47)  uses  only  the  first  k  iterates  of  the  generated  sequence, 
the  most  accurate  approximation  of  the  limit  for  a  convergent  sequence,  £fc+i,  is  not 
used.  What  would  be  helpful  is  to  modify  the  algorithm  so  that  the  most  current 
estimates  are  used.  Chandler  (1988)  suggests  the  following  model  for  the  linear 
equation  (4). 

Let  the  finite  sequence  5  =  {x0,xi, . . . ,  Xfc+i}  be  the  k  +  2  generated  vectors 
that  are  used  to  determine  the  coefficient  vector  c  of  (44).  As  shown  in  Chapter  II, 
Equation  (13),  there  exist  constants  a,'  and  i  =  1, . . . ,  m,  where  m  is  the  order  of 
A ,  such  that 

m 

Xn  =  ^2  di  Vi  q”,  n  =  0, 1, . . . ,  k  +  1. 

«=i 

Define  =  a^i+l  and  #  =  1 /$,*'  =  1, ...  ,m.  Then 
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x„  =  s  +  ^2  bi  Vi  Pi +1  n,  n  =  0, 1,  . . ,  k  +  1. 

t=i 

Then  there  exist  a  sequence  T  =  {x  '0,x  .  ,x'k+1}  such  that 

m 

x  '  =  i*+  Vj  p*,  r  =  0, . . . ,  k  +  1, 

t=0 

where  x  '  =  Xfc+i_r.  If  |  <  1  for  all  i,  then  |pj|  >  1.  Hence,  S' is  the  limit  of  S 
and  the  antilimit  of  T.  If  the  moduli  of  all  qC s  are  greater  than  unity,  then  |p;|  <  1 
for  all  i.  Therefore,  S'  is  the  antilimit  of  S  and  the  limit  of  T.  Otherwise,  S'  is  the 
antilimit  of  both  S  and  T. 

If  the  generated  sequence  S  diverges,  then  the  most  accurate  estimate  of  the 
antilimit  is  x0.  Therefore,  the  MPE  method  should  be  applied  to  the  sequence  S 
in  that  case  to  approximate  the  antilimit  so  that  (47)  will  be  the  sum  of  the  most 
accurate  estimate  plus  a  small  error.  If  S  converges,  then  the  most  accurate  estimate 
of  the  limit  is  Xfc+i;  hence,  sequence  T  should  be  used  in  this  case  to  approximate 
the  limit.  Even  though  the  above  theory  was  developed  for  the  linear  case,  it  can  be 
used  as  a  model  for  estimating  the  solution  of  a  nonlinear  problem.  A  comparison 
of  the  two  techniques  will  be  shown  in  Chapter  XI. 


CHAPTER  VIII 


THE  REDUCED  RANK  EXTRAPOLATION  METHOD 
Theory  for  the  Full  Rank  Extrapolation  Method 

Henrici  (1964)  set  forth  to  extend  Aitken’s  formula  for  systems  of  equations. 
His  goal  was  to  estimate  the  limit  of  a  sequence  of  m-dimensional  vectors.  His 
formula  contains  two  m  x  m  matrices,  of  which  one  involves  an  inverse.  For  large 
problems,  solving  a  large  linear  system  is  not  exactly  helpful.  However,  the  theo¬ 
retical  application  of  his  work  is  valid.  Meiina  (1977)  and  Eddy  (1979)  modified 
Henrici ’s  basic  formula  by  reducing  the  dimension  of  the  linear  system  to  a  value 
that  is  reasonable  for  computation.  Eddy  referred  to  his  method  as  the  Reduced 
Rank  Extrapolation  (RRE)  method. 

Before  their  methods  and  formulas  are  derived,  some  basic  definitions  are 
needed  which  will  be  used  throughout  this  chapter.  Some  of  the  definitions  have 
already  been  used  in  previous  chapters;  however,  they  are  mentioned  again  for 
completeness. 

Let  {®n}  be  an  m-dimensional  vector  sequence  generated  by  the  Equation  (3) 
such  that  S’  is  a  solution  of  (2).  Define  the  first  and  second  difference  vectors  as 
un  =  xn+i  —  xn  and  vn  =  un+1  —  un,  respectively.  The  following  m  x  k  (1  <  k  <  m) 
rectangular  matrices  are  very  valuable  in  the  development  of  the  theory.  Their 
columns  are  first  or  second  difference  vectors.  Define 

Un  =  [u„,u„+i,...,un+fc_1]  and 

=  [Vri5  ,  .  .  .  ,  Vn+j|,_i]  . 
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Also,  define  A  =  A(J)  to  be  the  Jacobian  matrix  of  the  function  G  taken  at  the 
solution  s.  The  Jacobian  matrix  is  defined  as 

9i(x  ) 

9m(x  )  _ 

Oij  is  the  (i,j)  component  of  the  matrix  A,  and  x  =  (®i,  ®2>  •  •  •  >  xm)T .  Henrici  (1964, 
p.  104)  showed  that  en+i  =  A(s)en  +  0(||e„||)2  where  en  =  xn  —  s  and  0(||en||2) 
denotes  a  quantity  bounded  by  C'||en||2,C  an  integer.  If  we  assume  that  this  error 
formula  is  exact  with  0(||en||)  =  0  for  finite  values  of  n,  then  xn+i  —  s  =  A(xn  —  s). 
Then  the  following  relationships  axe  satisfied: 


an  = 


dgijx  ) 
dxj 


*,j  =  1, . . .  ,m,  where  G(x  )  = 


un  =  Aun. i  =  •  •  •  =  Anu o,  vn  =  (A  -  I)un, 

Vn  =  Un+1-Un  =  (A-I)Un,  and  (48) 

Un+i  ~  AUn.  (49) 


If  k  —  m,  then  Un  and  Vn  are  square  matrices.  Assuming  that  Un  is  nonsin¬ 
gular,  Equation  (49)  gives 

A  =  Un+1U~\  (50) 

Since  A  is  the  Jacobian  matrix  of  G,  then 


®n+i  &  —  Axn  As. 


This  implies  that 


(A  f)>s  Axn  zn+i  -b  xn  xn  —  (A  J")xn  un. 
Therefore,  if  unity  is  not  an  eigenvalue  of  (A  -  /),  then  (A  -  I)-1  exists  and 

S=i„-(A-  I)"1 


(51) 
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Using  (50),  we  change  ( A  —  J)_1  in  the  following  fashion: 

(A  - 1)-'  =  (tUitr-1  -  /)-'  =  [(tr„+,  -  U„)v;'Yl 
=  [W-*]_1  =  u.v-K 

Substituting  into  (51),  the  extrapolation  formula  is 

s  =  xn  -  UnV~lun.  (52) 

Though  the  development  of  (52)  is  due  to  Henrici  (1964),  some  of  the  notation  used 
is  due  to  Smith,  Ford,  and  Sidi  (1987).  As  was  the  case  with  the  MPE  method, 
Equation  (52)  should  be  applied  to  the  sequence  S  =  {xn,x„+i, . . .  ,xn+fc+1}  if  the 
sequence  S  diverges.  If  S  is  a  convergent  sequence,  then  the  sequence 
T  =  {xn+fc+i ,  xn+fc, . . . ,  xn}  is  used  instead  of  S. 

Eddy  (1979)  derived  the  same  extrapolation  formula  as  follows: 

s  =  lim  xn  =  lim(x0  +  (xi  -  x0)  +  (x2  -  xi)  +  . . .) 

=  lim(x0  +  Uq  +  U\  +  . .  •) 

=  x0  +  lim  (I  +  A  +  A2  +  ..  .)u0 

=  X0  +  {I-  A)_1u0.  (53) 

Equation  (48)  then  yields 

s  =  x0  -  C/oF0-1tfo,  or 
s  —  xq  +  UqZ  and  0  =  tto  4*  Vo 

which  matches  Equation  (52)  for  n  =  0.  Therefore,  for  a  linear  system  and  with  no 
rounding  errors,  sis  computed  exactly.  For  a  nonlinear  problem,  the  limit  cannot  be 
obtained  for  a  finite  value  of  k,  but  Equation  (52)  is  used  as  a  model  for  estimating  s, 
the  solution  of  the  problem.  Since  the  extrapolated  vector  will  be  only  an  estimate 
of  3,  a  repeating  process  with  a  new  initial  vector  set  equal  to  the  extrapolated 
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vector  is  used  to  establish  a  semi-dynamic  procedure.  This  method  will  be  referred 
to  as  the  Full  Rank  Extrapolation  (FRE)  method. 

The  FRE  method  is  an  acceleration  technique  that  requires  m  +  1  iterations, 
where  m  is  the  dimension  of  the  vector  space  of  the  problem.  There  is  one  obvious 
problem  with  this  method:  if  m  is  large,  then  to  obtain  m  +  1  iterations  before 
we  can  even  apply  the  extrapolation  technique  defeats  the  purpose  of  accelerating. 
Though  Henrici  (1964)  indicated  that  the  technique  is  still  valid  for  values  much 
smaller  than  m,  Eddy  (1979)  proved  this  fact. 

Theory  for  the  Reduced  Rank  Extrapolation  Method 

Assume  that  we  choose  k  such  that  1  <  k  <  m.  Then  U  and  V  are  now  non¬ 
square  matrices  and  have  no  inverses,  so  Equation  (50)  does  not  hold.  Therefore, 
an  alternate  approach  for  establishing  the  basic  extrapolation  formula  is  needed. 
Let  the  exact  limit,  s,  in  (53)  be  replaced  by  the  extrapolated  value  s' .  Define 

(I  -  A)_1u0  =  U0Z.  (54) 

Using  (48)  and  (54),  we  have 


u0  =  {I  -  A)U0Z  =  -V0 Z 

so  that  the  extrapolated  vector  cam  be  expressed  by 

s  —  Xo  -f-  UqZ  and  0  =  uq  +  Vo Z. 
Solving  Equation  (55)  by  the  method  of  least  squares  gives 


Therefore, 


0  =  (Vq)uq  +  {{Vq)Vo)Z. 


(55) 


2  =  -(WJVcJ-'TOuo  =  -V0+u„, 


(56) 
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where  is  the  generalized  inverse  of  Vq. 

Substituting  (56)  into  (55)  and  generalizing  give  an  extrapolation  method  for 
k  <  m  and  any  starting  vector  £n: 

s  =  xn-UnV+un.  (57) 

As  with  the  MPE  and  FRE  methods,  the  sequence  used  in  (57)  may  be  the  generated 
sequence,  5,  or  S  in  reverse  order,  depending  upon  whether  the  iterated  sequence 
diverges  or  converges,  respectively.  Eddy  (1979)  called  this  method  the  Reduced 
Rank  Extrapolation  (RRE)  method.  Once  again,  it  is  usually  best  to  keep  k  less 
than  about  six  for  large  dimensional  problems.  From  this  point  on,  this  procedure 
will  be  referred  to  as  the  RRE  method,  regardless  of  the  value  of  k. 

REDUCED  RANK  EXTRAPOLATION  ALGORITHM  8.1: 

Given  the  iteration  equation  Xi+1  =  (7(£j),  the  initial  vector  x0,  and  the  positive 
integer  k. 

Step  1.  For  i  =  0, 1, . . . ,  k,  compute  ij+1  =  G(xt). 

Step  2.  For  i  =  0, 1, . . . ,  k  —  1,  compute  U{  =  xt+1  —  £  and  £;  =  xfi+i  —  Uj. 
Step  3.  Define  U  and  V  by  U  =  [u0,ui, . . .  ,Ufc_a]  and  V  =  [£0,  Vi,  •  •  •  ,vfc-i]- 
Step  4.  Compute  J*  =  xq  —  U Hu o,  where  H  =  V-1,  if  k  =  m,  the  dimension 
of  the  problem;  or  H  =  V+  =  ((V*)^)- 1 V*,  the  generalized  inverse 
of  V,  if  k  <  m. 

Step  5.  If  HG(i')  —  s\\  <  Tol,  then  stop;  otherwise,  set  x0  =  a  and  ii  =  G{s), 
generate  the  vectors  £2, . . .  ,£*+i  by  the  iteration  equation,  and  go 
to  step  2. 

The  computation  of  Step  4  of  Algorithm  8.1  involves  the  generalized  inverse 
of  V.  Eddy  (1979)  and  Smith,  Ford,  and  Sidi  (1987)  suggest  that  this  matrix  be 
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computed  in  the  algorithm.  Therefore,  the  matrix  ((F*)F)-1  must  be  determined, 
which  requires  costly  computer  time,  especially  for  larger  values  of  k.  I  suggest  am 
alternate  approach.  Let  B  =  (V*)V  and  y  =  ( V*)uo .  Then  UHxlq  in  Step  4  can  be 
rewritten  as 

UHu0  =  IT((V'*)V)“l(V’*)t^  =  UB~1y  =  Uz , 

where  Bz  ~  y.  The  vector  z  is  found  by  Gaussian  elimination.  Since  the  product 
(V*)V  will  be  a  symmetric  matrix,  the  amount  of  computer  time  is  reduced  even 
more. 


Numerical  Examples 

For  an  example,  consider  the  two-dimensional  problem  (Henrici,  1964)  of  find¬ 
ing  the  solution  of  the  system 

x  =  x2  +  y 2, 
y  =  x2  -  y2, 

near  the  point  (0.8, 0.4).  A  quick  check  will  show  that  the  solution  of  (58)  is 
(0.771845,0.419643)  to  six  decimal  places.  However  converting  (58)  into  its  it¬ 
eration  equations, 

sn+1  =  (®„)2  -1-  ( yn )2  and  y„+i  =  (sn)2  +  (yn)2, 

the  iterative  sequence  {zn}»  where  zn  =  ( xn,yn ),  diverges. 

Table  13  (page  67)  shows  infinity  norm  results  of  applying  the  RRE  method, 
k  —  2,  to  this  problem.  The  first  column  gives  the  norms  of  each  difference  vector, 
un.  The  norms  given  in  the  second  column  are  for  the  error  vectors,  e*n,  of  each 
extrapolated  vector.  Since  k  =  2,  there  must  be  k  +  1  =  3  iterations  before  the 
extrapolation  technique  can  be  applied.  Hence,  extrapolated  results,  column  two, 
are  obtained  only  once  every  three  iterations. 
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TABLE  13 


INFINITY  NORMS  OF  DIFFERENCE  AND 
ERROR  VECTORS  FOR  THE  RRE 
METHOD  FOR  PROBLEM  (58) 


n 

0 

1 

2 

3 

4 

5 

6 

7 

8 

9 


0.080000 
0.070400 
0.180224 
0.003916 
0.004952 
0.009096 
0.000023 
0  000047 
0.000079 
0.000000 


0.002280 


0.000023 


0.000000 


Table  13  shows  that  the  solution,  to  six  place  accuracy,  is  obtained  upon 
extrapolating  after  the  ninth  iteration.  However,  a  tenth  iteration  is  needed  to 
ensure  that  the  ninth  difference  vector  has  the  desired  precision  of  accuracy  for 
convergence.  In  addition,  column  two  measures  the  error  vector,  which  normally 
cannot  be  measured  since  the  answer  is  not  known. 

Figure  15  (page  68)  shows  a  graph  of  how  the  RRE  method  ( k  =  2)  compares 
with  the  MPE  method  ( k  =  2),  the  modified  vector  epsilon  method  (order  2),  and 
Aitken’s  semi-dynamic  method  (Jennings’  SDM)  on  Problem  (58).  The  graph  plots 
the  logarithm  (base  10)  of  the  infinity  norm  of  the  difference  vector 
(«„-!  =  xn  —  fn-i)  as  a  function  of  the  number  of  iterations.  The  results  show  that 
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the  RRE  and  MPE  methods  clearly  converge  faster  then  the  e  method  and  Aitken’s 
method.  In  fact,  the  graph  suggests  that  the  RRE  and  the  MPE  are  equivalent 
methods  since  the  results  obtained  from  these  two  methods  are  identical.  However, 
this  is  NOT  so.  In  Chapter  X,  it  will  be  shown  that  the  two  methods  are  very 
similar  (in  fact,  their  results  are  identical  for  some  problems,  as  is  the  case  for  this 
example),  but  they  are  not  equivalent  methods.  Probably,  the  most  important  fact 
shown  in  Figure  15  is  that  the  MPE  and  RRE  methods  both  produced  an  accurate 
fixed  point  solution  from  a  divergent  nonlinear  iteration  and  obtained  the  results 
very  rapidly. 


Figure  15.  Graph  Comparisons  for  Problem  (58) 


CHAPTER  IX 


ANDERSON’S  GENERALIZED  SECANT  ALGORITHMS 
Theoretical  Development  for  Secant  Methods 

Anderson  (1965)  was  motivated  by  the  inability  of  Aitken’s  and  Wynn’s  meth¬ 
ods  to  “feed  back”  into  the  process  in  a  fully  dynamic  manner  iterates  already  ob¬ 
tained.  He  expressed  his  feelings  by  stating  “The  Aitken  A2  process,  of  which  the 
e-algorithm  is  a  generalization,  is  considerably  less  effective  if  applied  statically  ... 
than  if  applied  dynamically  ...”  (Anderson,  1965,  p.  552).  His  referral  to  a  dynamic 
process  is  what  has  been  called  a  semi-dynamic  method  in  this  study.  He  desired  to 
find  a  fully  dynamic  procedure  which  would  accelerate  the  convergence  of  a  vector 
sequence.  This  process  is  similar  to  the  fully  dynamic  scalar  Aitken  algorithm  of 
Irons  and  Shrive  (1987)  discussed  in  Chapter  IV.  Anderson  developed  one  algorithm 
and  then  derived  variations  from  it.  He  obtained  the  first  algorithm  by  generalizing 
the  univariate  secant  method  geometrically;  see  Figure  16  (page  70).  Given  the 
equation  x  =  g(x)  and  the  scalars  xq  and  *i,  the  univariate  secant  method  is 

xn+ 1  =  *n  4"  -®(*n-l  —  *n)>  ^  =  1>  2,  .  .  .  , 

where 

q  _  _ 9{xn)  ~  xn _ 

[g(xn)  -  0(*„-l)]  -  [*»  -  *n-l]  ’ 

Generalizing  the  method  for  m-dimensions,  Wolfe  (1959)  saw  the  next  element  of 
the  sequence  as  being  the  solution  of  a  system  of  nonlinear  equations  of  m  secant 
hyperplanes  through  m  +  1  points.  However,  Anderson  considered  only  a  hyperline 


69 


70 


Figure  16.  Graph  of  One  Extrapolation  of  Secant  Method 


through  two  points.  It  .must  be  noted  that,  in  general,  a  hyperline  does  not  intersect 
the  subspace  defining  the  solution;  however,  the  point  chosen  is  the  point  that  is  in 
some  sense  “closest”  to  this  subspace. 

Now  for  the  development  of  Anderson’s  first  algorithm.  Given  the  vector 
sequence  {*„}  and  the  basic  iteration  equation  zn+i  =  G(xn),  Anderson  sought  two 
other  sequences  which  converge  to  the  same  limit  as  {:?„},  but  more  rapidly.  He 
first  defined  a  coupled  pair  of  iterative  sequences  {yn}  and  {z„}  by 

zn  =  G(yn). 

Also  define  the  residual  vector,  r*„,  and  the  inner  product,  (if,  if),  of  the  two  real 
m-dimensional  vectors  u  and  v  by 

m 

rn  =  zn  -  yn  and  ( u,v )  =  v<  w»> 

i  =  l 
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respectively,  where  the  weights  W{  are  positive  real  numbers  and  the  scalars  Ui  and 
V{  axe  the  components  of  u  and  v.  Define  u'n  and  v'n  for  the  generalized  univariate 
case 


u'n  =  Vn  +  Qn(yn-1  -yn),  and 

^  n  ~  Zn  "b  Qn(Zn— 1  ^n). 

Also  define  the  “linearized  residual”  Rn  by 


(59) 


Rn  =  0.5(v'n  -  U'n,v'n  -  U'r). 

If  R  =  ( x ,  x),  then  from  calculus  =  2  x  j .  Hence,  minimizing  Rn  with  respect 
to  the  parameter  Qn  yields 


'dv’n  du'n 


dRn 

dQn  \8Qn  8Q„ 

=  (fn_!  -rn,v'n-u'n)  =  0. 


Solving  for  Qn ,  we  have 


Hence, 


0  —  (rn-l  —  Tn,rn  +  Qn  [rn-l  ^*n]) 

=  (rn-l  -  rn,  rn )  +  Qn(rn- 1  —  rn>  rn- 1  ~  rn)- 

Qn  —  (b*n  —  *V»-l>rn)/(rn  —  rn-l,rn  ~  rn-l  )• 


(60) 


Define  the  extrapolated  vector  by 


&+i  =u'n  +  Bn(v'n-u'n),  Bn  >  0. 


(61) 


According  to  Anderson  the  choice  of  a  positive  Bn  prevents  yn+j  from  becoming 
trapped  in  the  subspace  spanned  by  the  previous  yn  iterates.  Usually  Bn  =  1  is  most 
appropriate;  however,  one  must  determine  the  optimum  value  for  Bn  empirically. 
Anderson  refers  to  this  algorithm  as  the  “extrapolation  algorithm.” 

Before  applying  Anderson’s  extrapolation  method,  the  first  two  terms  of  both 
{yn}  and  {£„}  are  required.  Therefore,  {jfn}  and  {£„}  usually  are  initiated  by 
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setting  yo  =  x0  and  computing  z0  =  G(yo)  =  Also  set  yi  =  z0  and  compute  zi 
by  zi  =  G(yi)  =  (7(xi)  =  The  extrapolation  technique  is  now  applied. 

For  the  special  case  where  m  =  1,  Qn  =  rn/(rn  —  fn-i).  Therefore,  for  Bn  =  0, 
the  next  iterate  is 


Vn 


2/n+l  2/n  A  /  -.  \  /-•  -»  \ 

(•^n  2/n)  (^n— 1  2/n—  l) 

Substituting  zn  =  G^yjj)  and  rearranging  give 

G{yn)  -  yn 


(j/n-l  2/n)> 


2/n+l  —  2/n 


Vn  —  yn- 1  —  G(yn)  +  G(yn-i) 


(iln  -  2/n- 1), 


which  is  the  univariate  secant  method.  Hence,  Anderson’s  method  is  consistent  for 
m  =  1. 

ANDERSON’S  EXTRAPOLATION  ALGORITHM  9.1 


Given  the  iteration  equation  xn+1  =  G(xn),  the  initial  vector  x0,  and  the  sequence 
. . .}. 

Step  1.  Define  y0  =  x0,  z0  =  yi  =  z\,  =  z0  —  y0?  and  set  n  —  1. 

Step  2.  Compute  zn  =  G(yn)  and  fn  =  zn  —  yn. 

Step  3.  Find  Qn,u'n,  and  v'n  by  (60)  and  (59). 

Step  4.  Compute  yn+1  =  u  'n  +  Bn{v  'n-u'n). 

Step  5.  If  ||2/n+i  —  yi* |j  <  Tol,  stop;  otherwise,  increase  n  by  one  and 
go  to  step  2. 

Anderson  developed  an  alternate  algorithm  he  referred  to  as  the  “relaxation 
algorithm.”  The  name  was  so  given  because  the  method  defines  a  relaxation  param¬ 
eter  dynamically.  Define 

u'n  =  zn  +  Qnrn, 

v'n  =  zn-i  +  <?nrn_i,  and 


Rn  =  0.5(u  —  u'n,v'n  —  £.*(,). 
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Minimizing  Rn  with  respect  to  Qn  yields 


Hence, 


dRn 

dQn 


(fn_!  -  rn,v'n-u'n )  =  0. 


0  —  (f*n-l  7*n>(^n-l  ?n)  +  Qn  [’’n-l  rn]) 

=  (j'n-l  —  rn,zn-l  ~  zn  +  Qn(rn- 1  —  rn,rn~ 1  —  rn)- 


So 


Qn  —  (^*n  J'n-li  2n_i)/(rn  T"n— lj^n  ^ n— 1)>  (62) 


Thus  define  yn+1  =  u'n.  It  should  be  noted  here  that  Anderson  has  a  typographical 
error  in  his  article.  His  equation  does  not  have  the  negative  sign. 


ANDERSON’S  RELAXATION  ALGORITHM  9.2. 


Given  the  iteration  equation  xn+i  =  G(xn),  the  initial  vector  z0,  and  the  sequence 

{Bo,  Bi,.. .}. 

Step  1.  Define  y0  =  x0,  zo  =  Vi  =  X\,fo  =  Zo  —  yo,  and  set  n  =  1. 

Step  2.  Compute  zn  =  G(yn)  and  rn  -  zn  -  yn. 

Step  3.  Determine  Qn  by  (62). 

Step  4.  Set  yn+1  =  u'n  =  zn  +  Qnrn. 

Step  5.  If  \\yn+i  —  yn||  <  Tol,  stop;  otherwise,  increase  n  by  one  and 
go  to  step  2. 

Anderson  discussed  two  particular  variants  of  the  first  algorithm.  The  first  is 
the  choice  of  the  metric  of  the  inner  product  for  which  Rn  is  defined.  The  second 
variant  is  for  higher  degree  methods.  The  higher  degree  methods  are  obtained  by 
minimizing  a  linearized  residual,  Rn,  over  subspaces  cf  higher  dimensions. 
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Define,  for  a  positive  integer  k , 

k 

U'n  =  yn  +  'EQniVn-j  -Vn), 

3= 1 
ft 

v'n  =  2  ~  *")>  and  (63) 

j=l 

J?n  =  0.5  (v'n-u'n,v'n-u'n). 

Minimizing  Rn  with  respect  to  Q\  yields,  for  i  =  1, 2, . . . ,  k, 

k 

rn- i,  rn  rn—j  )Qn  ~  (rn  rn—i,  rn)'  (®4) 

3=1 

Define  yn+i  as  in  Equation  (61). 

This  algorithm  is  dynamic,  coupled,  and  can  be  applied  after  finding  k  +  1 
iterations.  In  addition,  Anderson’s  method  can  “build  up”  the  degree  by  being 
applied  for  k  =  1,2,  etc.,  until  the  desired  value  of  k  is  reached.  Hence,  for  any 
positive  integer  k ,  Anderson’s  higher  degree  algorithm  is  a  fully  dynamic  technique 
which  can  be  applied  after  only  two  iterations.  As  previously  mentioned,  Anderson 
states  that  low-degree  cases,  limiting  k  to  less  than  six,  usually  work  best. 

ANDERSON’S  HIGHER-DEGREE  ALGORITHM  9.3. 

Given  the  iteration  equation  z„+i  =  G(xn),  the  initial  vector  x0,  the  positive  integer 
k,  and  the  sequence  {Bo,  B\,.. .}. 

Step  1.  For  n  =  0  to  k,  define  yn  =  xn,  zn  =  G(yn),  and 

Tn  ~  J/n* 

Step  2.  Solve  the  system  (64)  for  Q £. 

Step  3.  Determine  u'n  and  v'n  by  (63). 

Step  4.  Compute  yn+x  =  u  'n  +  Bn(v  'n-u'n). 

Step  5.  If  ||yn+i  —  yn||  <  Tol,  stop;  otherwise,  compute 

zn+ 1  =  G(y„+i),rn+1  =  zn+x  -  y„+i ,  n  =  n  +  1,  and  go  to  step  2. 
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Numerical  Examples 

Consider  the  linear  system  (4)  where  A  is  the  tridiagonal  matrix  whose 

superdiagonal  is  (1,0, 1/3, 0,0,0), 

diagonal  is  (1/2, 1/2, 1,-1, -1/6, 1/3, 1/3, 1/3), 

(65) 

subdiagonal  is  (0,0, —1,0, 3,0),  and 

b  is  the  constant  vector  (—1/2, 1/2, —1/3, 13/6,2/3, —7/3, 2/3)r. 

Figure  17  compares  the  convergence  of  Anderson’s  method  for  2  <  k  <  5  and  the 
relaxation  method  with  the  basic  iteration.  Convergence,  the  infinity  norm  of  the 
difference  vector  less  than  10~1S,  with  no  acceleration  is  obtained  in  57  iterations. 
The  best  results  are  obtained  for  Anderson’s  method  with  k  =  4  and  5.  This  is 
not  by  coincidence  since  this  problem  was  designed  to  have  a  matrix  of  order  seven 
with  a  monic  minimal  polynomial  of  degree  4: 


Figure  17.  Graph  Comparison  of  Anderson’s 
Methods  for  Problem  (65) 
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Figure  18.  Comparison  of  Results  for  Anderson’s  Methods, 
the  MPE  Method,  and  the  RRE  Method  for 
Problem  (28) 


p(y)  =  (y-  1/2 )2(y  -  1/3)2, 

Hence,  the  best  results  are  obtained  when  k  is  chosen  as  the  degree  of  the  minimal 
polynomial  of  A,  as  was  the  case  for  the  MPE  and  RRE  methods  for  linear  problems. 
Another  point  is  that  if  k  is  greater  than  the  degree  of  the  minimal  polynomial,  the 
sequence  will  still  converge;  however,  the  convergence  rate  will  not  improve.  Once 
again,  we  see  that  an  acceleration  method  has  determined  a  fixed  point  solution 
much  faster  than  the  basic  iteration. 

Recall  the  familiar  examples,  Problems  (28)  and  (58).  In  Figures  18  and  19 
(Figure  19  on  page  77),  the  results  of  these  two  problems  for  Anderson’s  three 
methods:  k  =  1,2,  and  the  relaxation  method,  are  compared  with  the  results  ob¬ 
tained  by  the  RRE  and  MPE  methods,  k  —  2.  Clearly,  Anderson’s  method  with 
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Figure  19.  Comparison  of  Results  for  Anderson’s  Methods, 
the  MPE  Method,  and  the  RRE  Method  for 
Problem  (58) 


k  =  2  gives  the  best  results  for  these  examples.  This  leads  to  a  few  questions.  Are 
these  results  consistent  with  results  of  other  problems?  Are  there  ways  of  comparing 
two  or  more  of  these  techniques  from  the  theoretical  view  point?  The  latter  question 
leads  into  the  next  chapter  where  the  first  extrapolation  method  of  Anderson,  the 
MPE  method,  and  the  RRE  method  are  compared  in  theory.  The  former  question 
is  saved  for  Chapter  XI. 

A  final  comment  concerning  Anderson’s  method  needs  to  be  made.  Even 
though  he  wrote  his  article  in  1965,  it  has  been  widely  ignored.  There  were  no 
references  found  to  Anderson’s  article  or  to  Anderson’s  method  in  the  research  for 
this  thesis.  The  majority  of  the  articles  or  papers  written  on  the  subject  area  concern 
the  theoretical  and/or  numerical  comparison  of  the  vector  Aitken,  the  vector  c,  the 
RRE,  and/or  the  MPE  methods.  One  may  speculate  that  the  title  of  his  article, 
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‘Iterative  Procedures  for  Nonlinear  Integral  Equations,”  may  have  something  to  do 
with  the  problem  since  it  makes  no  reference  to  the  subject  of  acceleration  methods. 
For  whatever  the  reason,  the  article  and  the  acceleration  method  have  received  little 


attention. 


CHAPTER  X 


THEORETICAL  COMPARISONS 

Determinant  Form  for  the  MPE  Method 

Now  that  all  the  acceleration  techniques  have  been  given,  they  can  be  com¬ 
pared  theoretically.  Since  the  MPE,  the  RRE,  and  Anderson’s  methods  all  rely  on 
solving  a  system  of  linear  equations,  it  seems  logical  to  start  with  them.  The  theory 
presented  in  this  chapter  for  the  MPE  and  the  RRE  methods  was  derived  by  Sidi 
(1986)  and  Sidi,  Ford,  and  Smith  (1986).  The  theory  for  Anderson’s  method  and 
the  examples,  to  the  best  of  our  knowledge,  have  not  been  published. 

Define  the  inner  product  of  two  sequence  terms  a+  and  a,j  to  be  a,i-7  =  (a^,  dj), 
i,j  >  0.  Also  define  the  matrix  [t0,  . . . ,  tfc]  by 


to 

tl 

tfc 

0-0,0 

<*0,1 

Oo.fc 

<*1,0 

<*1,1 

<*1,* 

(66) 

Ofe-1,0 

<*fc-l,l 

<*fc-l,Jt 

Denote  the  determinant  of  the  matrix  (66)  by  D  [t0,  •  •  .,<*],  and  denote  by  JV,-  the 
cofactor  of  tj  in  22  [to, . . . ,  t&],  i  =  0,  ...,&.  If  the  elements  of  the  first  row  of 
D  [to,  •  • .  are  vectors,  then  the  determinant  is  to  be  interpreted  as 

k 

D  [f0,...,ffc]  =  ££  Ni.  (67) 

»=o 
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It  was  shown  in  Chapter  VII  that  given  the  m-dimensional  vector  sequence 
{zn},  its  limit  S',  and  some  positive  number  k  <  m;  a  can  be  estimated  by  use  of 


the  MPE  method  by 


k  /  fe  ^ 

"  =  2  where  U  =  c;/  ^  c3-  , 

1=0  \i=o  J 


provided  Ej=oci  ±  The  vector  c  =  (c0,  cu . . . ,  ck-i)T  solves  the  system  of  equa¬ 
tions  Uc  =  —Uk,  where  U  and  uk  are  found  by  (45)  and  (38)  and  c*  =  1.  The  vector 
c  that  satisfies  this  system  satisfies  the  normal  equations 

J2(ui,Uj)cj  =  ~(ui,uk),  0  <  i  <  k  -  1, 

j=o 

Consequently,  the  vector  T=  (l0, ...  ,4)  of  Equation  (f)8)  satisfies  the  equations 

ZU‘i  =  h  “d  (69) 

=  0,  o<»<*-i, 

provided  these  equations  have  a  solution.  The  matrix  of  Equations  (69)  is  (66) 
where  U  =  1  and  <nj  =  Assuming  that  its  determinant  is  nonzero,  then 

Cramer’s  rule  can  be  used  to  write  the  solution  of  (69)  as 


,  =  Hi  =  o  <  7  <  Jb 

J  E,-=0 Ni  °'J_ 


From  (67)  and  (68)  we  find 

h  k  n- 

=  . 

_  Sj=o  ^ j  _  -P  [go>  •  •  •  >  (70) 

P[l,...,l]  "  2> [1 . 1]  * 

Determinant  Form  for  the  RRE  Method 

Now  consider  the  RRE  method.  From  Equation  (57)  and  letting  n  =  0,  the 
first  extrapolated  vector  found  by  the  RRE  method  is 


s'  =  z0-  U0V0+uo, 
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where  V+  is  the  generalized  inverse  of  V .  (If  k  =  m,  then  V+  is  the  inverse  of  V .) 
Define  q  =  (?o>  9i>  •  •  •  >  9fc-i)r  to  be  the  vector  which  satisfies  the  system  of  equations 


Vq  =  -u0. 


(72) 


Then  (71)  can  be  rewritten  as 


k- 1 


s'  =  xo  +  Uoq  =  z0  +  ^  q.:  tf,. 


(73) 


«=o 


As  was  the  case  for  the  MPE  method,  the  q  that  satisfies  (72)  will  also  satisfy  the 
normal  equations 

k- i 


Y^(vi,Vj)qi  =  -(t7;,u0),  0<i<k-  1. 
j=o 


(74) 


Substituting  vj  =  Uj+ 1  —  Uj  for  the  second  component  on  the  left-hand  side  of  (74) 
and  rearranging  the  equation,  we  have 


fe-i 


0  =  (tTi,tZ o)  +  -  Uj)qj 


3= o 


=  (Vi,Uo)  +  (Vi,Ul  -  Uo)qo  +  -  Uj)Uj 


3=1 


k- 1 


=  (vi,Uo)  +  (Vi,ui)q0  -  (vi,u0)qo  +  ^  [(»«>“>+ i)?i  ~  (v;,Uj)?j] 

3=1 

=  (vi,u o)(l  -  9o)  +  {vi,Ui){qo  -  qi)  +  ( Vi,u2){qi  -  q2 ) 

+  •••  +  (vi,Uk-i)(qk-2  -  qk- 1)  +  (vi,Uk)qk- 1 

fc-1 

=  (vi,tfo)(l  -  9o)  +  {vi,Uk)qk-i  +  ^2{vi,Uj){qj-i  -  fc),  0  <  i  <  k  -  1. 

3=1 


Define 


It  is  easy  to  see  that 


h 3  =  1  —  qo,lk  =  qk-i,  and 

1 3  =  qj-1  ~  ft,  1  <><*-!. 


Et  =  ‘- 

j=o 


(75) 


(76) 
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Therefore,  Equations  (75)  and  (76)  establish  a  one-to-one  correspondence  between 
the  ft’s  and  7,’s.  Hence,  the  system  of  linear  equations  (74)  for  the  ft’s  is  equivalent 
to  the  linear  system 

E?=0^  =  l>  and  (?7) 

=  o  <  *  <  —  l, 

for  the  Ij's.  Substituting  Uj  =  Xj+ 1  —  £j  on  the  right-hand  side  of  (73),  rearranging, 
and  applying  the  one-to-one  correspondence  (75),  (73)  can  be  written  in  the  form 
of  (68).  The  matrix  of  the  system  (77),  once  again,  is  (66),  where  a,j  =  (ft,u,). 
Assuming  that  the  determinant  of  the  matrix  is  nonzero,  Cramer’s  rule  can  be  used 
to  solve  the  system.  The  result  will  be  (70).  It  is  important  to  notice  that  even 
though  the  extrapolated  results  for  both  the  MPE  and  RRE  methods  have  the  same 
determinant  form,  the  results  axe  not  necessarily  the  same.  This  is  because  the  Oij 
elements  of  the  matrices  are  different.  For  the  MPE,  they  are  (ui, Uj)‘,  whereas,  for 
the  RRE,  they  axe  (ft,  Uj). 


Determinant  Form  for  Anderson’s  Methods 

Lastly  we  consider  Anderson’s  higher  degree  method.  Assume  that  the  first 
extrapolation  is  not  performed  until  k  +  1  iterations  have  been  obtained.  Hence, 
the  first  k  4-  1  iterations  will  generate  the  same  sequence,  x0, . . .  ,Xfc+i,  as  the  MPE 
and  RRE  methods  do.  Also  assume  that  B„  =  1  for  all  n.  Therefore,  using  (61) 
and  (63)  the  extrapolation  formula  becomes 

y„+l  =  v’„  =  zn  +  -£  QHz„.,  -  ?„).  (78) 

;=t 

As  was  shown  in  Chapter  IX,  to  start  Anderson’s  method  one  sets  yo  =  xo,z0  = 
G{y0)  =  xuui  =  *o,  etc.,  until  the  necessary  number  of  sequence  elements  are 
obtained  for  extrapolation.  Therefore,  the  two  sequences  yo, ...  ,yn  and  xq,  . . . ,  xn 
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are  the  same  sequences,  as  are  the  two  sequences  zq,  . . . ,  zn  and  xx, . . . ,  £n+i-  This 
information  will  be  very  valuable  later. 

The  Q:n' s  are  found  by  solving  the  linear  system  of  equations  (64).  Since 
fn  =  zn—yn  =  xn+i  —  xn  =  un,  then  (64)  is  equivalent  to 

k 

^(•un  -  un-i,un  -  un-j)qj  =  (un  -  un_„un),  1  <  i  <  k,  or 

j=i 

k 

^  )(^n  ^n— »— 1 ,  Un  'Un_j)^lj  =  ( Un  ,  tin) ,  0  <  l  ^  fc  1,  Or 

i=i 

fc-i 

^(un  -  =  (ftn  ~  £«),  0  <  Z  <  fc  -  1, 

j=0 

where  qj  =  The  system  can  be  rearranged  as  follows: 

0  =  -(U„  -  ■Un.i-ijZTn  -  Un-i)^  +  (u„  -  Un-i-'.,Un)  - 

k- 1 

^  ^ (zZn  ,  1Zn  ZZn_j_i  )9j  +  l 

J=1 

=  (zZn  ,  Un'jqi  +  (tin  Hn— t— 1 ,  1Xn)  +  (un  ZXn_j_i ,  Un_i  )qi 

fc-1  fe-1 

£(««  _  izn_i.1,Un)9i+1  +  -  ttn-i-l,«n-j-l)?i+l 

3-1  3=1 

k 

—  ( un  zzn_j_i ,  zin)(l  qi  q2  ■  •  •  9fc)  4~  ^  '(zzn  tzn_j_i ,  un_ j )qj . 

3  =  1 

Set 

Jo  =  (1  -  $1  -  ?a - and  1,-  =  1  <j<k. 

It  is  easy  to  verify  that  /,-  =  1.  Once  again,  we  have  established  a  one-to-otte 
correspondence  between  the  qS s,  0  <  t  <  k,  and  the  Ij’s,  0  <  j  <  k.  Hence,  the 
linear  system  of  equations  (64)  for  the  ©’ s  is  equivalent  to  the  linear  system 
Y*  l  •  =  1 

Sj=o(^n  l^n—i—  1  >  Ztn_j  )/y  =  0,  0  <  i  <  £  1, 

for  the  l/s.  As  with  the  other  two  methods,  the  matrix  for  this  system  of  equations 
is  (66)  where 

aiJ  =  (rn  —  rn-i-l  >  rn-j) 


83 


are  the  same  sequences,  as  are  the  two  sequences  z0, . . . ,  zn  and  £i, . . . ,  xn+i .  This 
information  will  be  very  valuable  later. 

The  Qi’s  are  found  by  solving  the  linear  system  of  equations  (64).  Since 
fn  =  zn  —  yn=  xn+i  —  xn  =  Uni  then  (64)  is  equivalent  to 

k 

£(un  ~  Un-i,Un  -  Un_i)9j  =  (un  ~  Un-i,U„),  1  <  *  <  &,  OT 

i=i 

k 

-  Un-i-uUn  ~  un-j)qj  =  (un  -  un-i-i,un)i  0  <  i  <  k  -  1,  or 

k- i 

^2(un  ~  un-i-\iUn  -  un-j-i)qj+i  =  (un  -  un.<-i,tt„),  o  <  i  <  k  -  1, 
i= o 

where  qj  =  Qn.  The  system  can  be  rearranged  as  follows: 

0  =  ~{un  -  Un-i-l,Un  -  Un-\)q\  +  (un  -  Un_i_i,Un)  - 
k-1 

XX**"  ~  Un-i-l,Un  ~  Un-j-l)qj+l 
j= 1 

=  ~(un  ~  Un-i-l  1  Un)qi  +  (u„  -  ,  Un)  +  (u„  -  Un-i- 1 ,  Un-1  )?1  ~ 

*-l  fc-1 

XX**"  -  Un-i-l,  «n)9j+l  +  XX**"  “  Un-i-l  i^n-j-l  )qj+l 
j= 1  j=l 

k 

=  (u„  -  tf„_i_1,un)(l  -  qi-  q2 - qk)  +  XI(**"  -  u„-,-i,  un-j)9i- 


/0  =  (1  -  9i  -  92 - 9fc)  and  l,  =  9,-,  1  <  j  <  k. 

It  is  easy  to  verify  that  £j=1  =  1.  Once  again,  we  have  established  a  one-to-otie 

correspondence  between  the  q^s,  0  <  i  <  k,  and  the  lj' s,  0  <  j  <  k.  Hence,  the 
linear  system  of  equations  (64)  for  the  9i’s  is  equivalent  to  the  linear  system 

£}=*«,  =  1, 

E‘=o(»n  -  “n-i-I.Un-jJO  =  0  <  1  <  *  -  1, 

for  the  lj' s.  As  with  the  other  two  methods,  the  matrix  for  this  system  of  equations 
is  (66)  where 

ai,j  =  (rn  —  rn-i-l ,  rn-j) 
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—  ( un  un-i-hun-j)- 


(79) 


Therefore,  the  solution  is 


h  = 


D[l, 


0  <  j  <  k. 


(80) 


Rearranging  (78)  and  using  the  one-to-one  correspondence,  we  have 

k 

2/n+l  =  v  n  ZZ  Zn  (  qj(zn-j  —  zn ) 

;= 1 
k 

=  ®n+l  +  2  4#(*n+l-i  -  *«+l) 
j=l 

=  (1  —  ?1  —  •  *  •  —  ?fc)®n+ 1  +  9l*n  +  •  •  •  +  qkXn+l-k 


=  Jo*n+l  +  ^l*n  +  h^n-l  +  '  ’  *  +  ?fc®n+l-fc 
k 

=  y~!  h  gn+l-j- 

J=  0 


Therefore,  from  (67)  and  (80),  and  letting  ra  =  fc,  we  have 

_  ^  iVjifc+i-j 

2/fc+l  v  *fc+l-i  r)  [1  1 1 

j=0  J=0  L1?  •  •  •  ?  AJ 

_  D  i  •  •  •  >  gi] 


(81) 


Anderson’s  extrapolated  vector  with  Bn  =  1  is  quite  different  from  the  ex¬ 
trapolated  vector  for  the  MPE  and  RRE  methods.  First  the  matrix  elements  aij 
are,  as  before,  different;  but,  in  addition,  we  see  that  the  first  row  of  the  matrix 
is  different.  Instead  of  the  sequence  elements  i0,  Xj, . . . ,  x*;  the  row  contains  the 
elements  xj{, . . . ,  x^» 

Consider  Anderson’s  higher  degree  method  with  one  change:  let  Bn  =  0  for  all 
n.  Anderson  said  that  this  value  should  never  be  used,  it  is  used  here  for  theoretical 
purposes  only.  Then  the  extrapolated  vector  will  be 

k 

Vn+l  =  u'n  =  yn  +  J2  Q3n{yn+1  -  jf«). 
j=l 


(82) 
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Since  all  other  factors  remain  the  same  as  for  the  case  where  Bn  =  1,  (82)  can  be 
rewritten  as 

k 

y„+i  =  u'n  =  yn  +  Yl  qj(yn-j  -  Vn) 

3=0 


=  xn  +  -  Xn)  4-  92(®n-2  -  In)  4 - h  $n(Zn-fc  -  In) 

=  (1  —  ?1  - - ?fc)*n  +  9l*n-l  +  ?2®n-2  + - H 

=  +  ^l*n-l  +  hXn-2  + - h  ?fc*n-fe 

=  (83) 

3=0 


As  was  done  in  deriving  (81),  (83)  is  used  to  obtain  the  extrapolation  formula 

_  D  [®fc,  •  •  ■  ,  2?o] 

Vk+1  -  D  [1, . . . ,  1]  ’ 

where  the  elements  a+j  of  Matrix  (66)  are  defined  by  (u/,  —  Uk-i~i,Uk-j)-  Hence, 
for  Anderson’s  two  cases,  the  extrapolated  vectors  have  the  same  determinant  form 
with  the  exception  of  the  first  row  of  the  numerator  matrices. 

Anderson  vs  RRE  Comparison 

Let  R  be  the  Matrix  (66)  defined  for  the  RRE  case,  and  let  A  be  the  matrix 
defined  for  Anderson’s  case  with  Bn  =  0  .  By  interchanging  the  columns  of  A,  A  can 
be  rewritten  as  (xo>  •  •  •  >  where  the  new  aij  elements  axe  ( u/,  —  Uk-i-x,Uj).  Since 
the  interchanging  of  columns  will  be  identical  for  both  the  numerator  and  denom¬ 
inator  matrices,  the  sign  change,  if  any,  will  cancel  out.  Hence,  the  extrapolated 
vector  has  not  changed. 

Since  the  two  matrices,  R  and  A ,  have  identical  second  components  for  the 
a.itj  elements,  let  us  consider  only  the  first  components  which  are  determined  by  the 
variable  i,  the  row  variable.  By  interchanging  rows,  we  can  rewrite  A  in  such  a  way 
that  the  jth  column  of  A  will  have  the  form 
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(uk  -  Uo ,Uj) 


(uk  ~  Uk-l,Uj). 

Once  again,  the  extrapolated  vector  will  remain  the  same  since  the  sign  change,  if 
appropriate,  of  the  numerator  and  the  denominator  will  cancel  out.  In  order,  for 
i  =  2  to  k  —  1,  set  the  ith  row  equal  to  the  ith  row  minus  the  (i  +  l)th  row.  The 
resulting  jth  column  of  matrix  A  will  be 

Xj 

(Hi  -  Uo,  Uj) 

(U2  -  Ul,Uj) 


(uk  -  Uk-l, Uj). 

This  column  is  identical  to  the  jth  column  of  R.  Therefore,  the  extrapolated  vector 
yit+i  is  the  same  vector  for  the  RRE  method  and  this  special  case  of  Anderson’s 
higher  degree  method,  Bn  =  0. 

Now  consider  Equation  (81)  for  Anderson’s  method  with  Bn  =  1  applied  to 
the  linear  equation  (4).  Let  a.j  =  Nj/D  [1, . . . ,  1] ,  j  =  0, . . . ,  k.  Then 

h  k  /  h 

Vh+1  =  53  ai  Xh+i-j  =  53  a3  A£h-j  =  A  53  ai 

j=0  j= 0  \j=o 

D  {x^ , . . . ,  Xo] 

D[  1 . 1]  ’ 

which  is  one  iteration  of  equation  (84).  Therefore,  it  is  seen  that,  for  the  linear 
case,  the  extrapolated  vector  for  Anderson’s  method  with  the  B's  equal  to  unity  is 
identical  to  one  iteration  of  the  extrapolated  vector  for  Anderson’s  method  with  the 
B's  equal  to  zero;  hence,  it  is  equivalent  also  to  one  iteration  of  the  extrapolated 
vector  obtained  by  the  RRE  method.  An  example  will  follow  to  illustrate  this  fact. 
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Numerical  Examples 


Smith,  Ford,  and  Sidi  (1987)  suggested  using  the  Gauss-Seidel  iteration  scheme 
applied  to  the  linear  equation  (Wynn,  1962,  Eq.  (14))  F(x)  =  Ax  —  6  =  0,  where 


.  2 

1 

3 

4 

10 

A  = 

1 

-3 

1 

5 

and  6  = 

4 

3 

1 

6 

-2 

8 

4 

5 

-2 

-1 

6 

(85) 


To  define  the  Gauss-Seidel  scheme,  let  cnj ,  t,  j  =  1, . . .  ,m,  where  m  is  the  order  of 
A,  be  the  (i,  j)  component  of  A.  Define  the  matrices  L  and  UP  by 


0  0  0  0 

a2i  0  •  •  •  0  0 

•  •  •  •  • 

Orol  OmJ  '  ’  *  ®m,ra- 1  0 


(86) 


UP 


0  di2  .  Ulm 

0  0  .  a2m 


0  0 


Om—  l,m 

0 


(87) 


Also  define  D  to  be  the  diagonal  matrix  [011,022, . . .  Then  the  Gauss-Seidel 

iteration  is  defined  by 

(D  +  L)xn+l  =  -(UP)xn  +  b.  (88) 


When  this  scheme  is  applied  to  Problem  (85),  the  result  is  a  divergent  sequence; 
however,  all  three  acceleration  techniques  obtain  the  solution,  (1,1, 1,1).  If  we  let 
k  =  2,  then  the  first  extrapolated  vectors,  to  five  place  accuracy,  for  the  RRE 
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method,  r,  and  Anderson’s  method,  a,  with  Bn  =  1,  are 

f  =  (-0.17247, 1.09243, 0.67697, 1.39913)r  and 
a  =  (0.64007, 1.43756, 1.24008, 1.26795)r,  respectively. 

Applying  an  iteration  to  r  shows  that  a  =  F(r),  as  was  shown  and  previously 
discussed.  This  fact  in  itself  would  lead  one  to  believe  that  for  linear  problems, 
Anderson’s  method  has  a  major  advantage  over  the  RRE  method.  Now  consider 
another  example,  a  nonlinear  one. 

In  Chapter  VIII,  Problem  (58)  was  used  to  compare  the  RRE  method  with 
previously  derived  techniques.  As  noted  in  the  last  paragraph  of  Chapter  VIII,  the 
RRE  method  and  the  MPE  method  seemed  to  be  equivalent  procedures  since  they 
produced  identical  results  for  k  =  2.  One  may  wonder  how  this  can  be  true  since  the 
theory  proven  in  this  chapter  shows  that  they  are  not  equivalent.  Let  us  examine 
this  problem  more  closely.  The  second  and  third  rows  of  matrix  (66)  for  the  MPE, 
RRE,  and  Anderson’s  methods  are  given  in  Table  14  (page  89).  Hence,  matrix 
(66)  is  different  for  each  method.  However,  when  the  extrapolated  vector  for  each 
method  is  computed,  we  find  that  all  four  methods  (MPE,  RRE,  and  Anderson’s 
method  with  Bn  =  0  and  Bn  =  1)  give  identical  results: 

(0.774124, 0.419430)r  =  1.080567xo  +  0.286985xa  -  0.367552x2 
for  the  RRE,  MPE,  and  Anderson’s  (i?n  =  0)  methods  and 

(0.774124, 0.419430)r  =  1.080567;?!  +  0.286985x2  -  0.367552x3 
for  Anderson’s  method  ( Bn  =  1),  where 


*o  =  (0.8, 0.4)r,  X!  =  (0.8,0.48)r, 

x2  =  (0.8704, 0.4096)r,  and  x3  =  (0.925368, 0.589824)r. 
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TABLE  14 

SECOND  AND  THIRD  ROWS  OF  MATRIX  (66)  FOR  THE  MPE, 
RRE,  AND  ANDERSON’S  METHOD  FOR  PROBLEM  (58) 


MPE  RRE 


0.00640 

-0.00563 

0.01442 

0.00802 

-0.00319 

0.02108 

-0.00563 

0.00991 

-0.00882 

0.02005 

-0.01873 

0.04432 

ANDERSON 

-0.01203 

0.01554 

-0.02324 

0.02005 

-0.01873 

0.04432 

Equivalent  results  were  expected  for  Anderson’s  method  with  Bn  =  0  and  the 
RRE  method;  but,  why  did  all  four  methods  obtain  the  same  results?  Consider,  for 
the  moment,  only  the  MPE  and  RRE  methods,  and  the  difference  between  their 
aitj  elements  of  matrix  (66).  Substituting  un+1  —  un  for  vn  and  carrying  out  the 
determinant  calculations,  the  resulting  coefficients  for  x0,xi,  and  x2,  respectively, 
in  terms  of  inner  products  with  notation  M0112  =  (u0,  u2),  are 

M1122  -  M0122  +  M0112  -  M 1212  +  M0212  -  M0211, 

M0212  -  M 0202  +  M0102  -  M 0122  +  M0022  -  M0012,  and 

M0112  -  AT0012  +  AfOOll  -  M0211  +  M01G2  -  JkfOlOl 

for  the  RRE  method  and 

M0112  -  M0211,  MO  102  -  M0012,  and  M0011  -  M0101 

for  the  MPE  method.  Comparing  the  different  coefficients,  there  is  a  difference  of 
M1122  -  M0122  -  M1212  +  M0212, 

M0212  -  M0202  -  M0122  +  M0C22,  and  (89) 

M0112  -  M0012  -  M0211  +  M0102 
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for  f0,xi,  and  x2,  respectively.  In  addition,  the  difference  of  D  [1, . . . ,  1]  for  the  two 
extrapolations  is 

M1122  -  2M0122  +  M0112  -  M1212  +  M0212  -  M0211+ 

(90) 

M0212  -  M0202  +  M0102  4  M0022  -  MO 012. 

Define  wQ,wl,w2  to  be  the  quotient  of  the  expressions  in  (89)  divided  by  (90)  for 
x0,  xu  and  x2,  respectively.  The  results  are 

wO  =  1.080567,  wl  =  0.286985,  and  w2  =  -0.367552. 

These  values  are  the  same  values  as  their  corresponding  coefficients. 

Consider  the  coefficient  of  one  term  of  D  [ x0,Xi ,  x2]  for  the  MPE  method.  Let 
NM  represent  this  value  and  let  DM  represent  the  value  D  [1,1,1].  In  addition, 
define  DN  as  the  difference  (89)  for  the  chosen  coefficient  and  DD  as  the  dif¬ 
ference  (90).  As  was  shown  in  the  previous  paragraph,  the  quotients  NM/DM 
and  DN/DD  are  equal.  Setting  this  value  to  r,  we  have  DN  =  DD(r)  and 
NM  =  DM(r).  Also, 

NM  +  DN  =  DM{r)  +  DD(r)  =  (DM  +  DD)r. 

Therefore,  (NM  +  DN)/(DM  +  DD)  =  r.  Since  this  is  true  for  each  coefficient, 
we  see  that  the  linear  combination  of  xVs  is  the  same  for  both  methods. 

Similar  results  can  be  shown  for  Anderson’s  method  with  Bn  =  1.  Though 
the  results  showed  equality  for  this  case,  it  is  easy  to  see  that  the  two  methods  are 
not  equivalent.  To  show  this  fact  by  example,  consider  the  problem 

x  =  x2  +  y2  —  z2, 
y  =  x2  -  y2  +  z2,  and 
z  =  —x2  ■+■  y2  4-  z2 . 


One  extrapolation  of  this  problem  gives 
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(0.175810,-0.786335,-0.978136)  for  the  MPE  method, 
(0.071951,0.020600,-0.039942)  for  Anderson’s,  Bn  =  1, 

(0.236638,0.039550,  — 0.063152)  for  the  RRE  method,  and 
(0.236638,0.039550,-0.063152)  for  Anderson’s,  Bn  =  0. 

Hence,  they  are  all  different  with  the  exception  of  the  RRE  and  Anderson’s  with 
Bn  =  0  methods,  as  expected.  The  above  results  indicate  that  the  MPE,  the 
RRE,  and  Anderson’s  ( Bn  =  0)  methods  are  very  similar  for  one  extrapolation. 
However,  the  first  extrapolated  vector  for  Anderson's  method  ( Bn  =  1),  for  the 
linear  case,  is  one  iteration  better  than  the  RRE  method.  If  the  fact  that  Anderson’s 
method  is  fully  dynamic  is  also  considered,  then  one  might  assume  that  Anderson’s 
method  ( Bn  >  0)  will  accelerate  a  sequence  to  its  correct  limit  faster  than  the  other 
techniques.  In  the  next  chapter,  we  will  check  the  validity  of  this  assumption  by 
testing  Anderson’s  method  and  the  other  methods  on  several  types  of  test  problems. 


CHAPTER  XI 


NUMERICAL  TEST  PROBLEMS 

In  this  chapter  we  compare  the  acceleration  methods  numerically.  Test  prob¬ 
lems  will  include  linear  and  nonlinear  problems  with  dimensions  varying  from  4  up 
to  8000.  The  problems  were  tested  using  FORTRAN  coded  programs.  All  problems 
except  Examples  9  and  10  were  computed  on  a  Kaypro  286  PC  with  double  pre¬ 
cision,  giving  a  relative  precision  of  1  x  10-19.  Examples  9  and  10  were  computed 
on  a  IBM  3081K  (VS  FORTRAN)  with  double  precision,  giving  a  relative  precision 
of  2  x  10-16.  Dr.  John  P.  Chandler,  Oklahoma  State  University,  developed  the 
main  software  for  the  semi-dynamic  Vector  Aitken  (VA)  method,  Wynn’s  original 
Vector  £  (Ve)  method,  the  Modified  Vector  e  (MVe)  method,  the  MPE  method, 
and  Anderson’s  (AND)  method.  I  wrote  the  program  for  the  RRE  method  and 
made  some  modifications  to  Chandler’s  Ve,  MPE,  and  AND  programs. 

The  results  will  be  presented  in  figures  which  show  graphs  of  the  logarithm 
(base  10)  of  the  infinity  norm  of  the  difference  vector,  un_i  =  xn— x„-i,  as  a  function 
of  the  number  of  iterations.  Exceptions  are  Examples  9  and  10.  Example  9  is  the 
graph  of  the  infinity  norm  ot  the  error  vector,  en  =  $  —  xn,  where  s  is  the  solution. 
Example  10  plots  the  Euclidean  norm  of  the  difference  vector.  The  notation  will  be 
logio  ||«n— 1||-  The  stopping  criterion  on  all  programs  is  ||u„||  <  1  x  10-15.  Therefore, 
this  value,  denoted  by  CT5,  will  be  considered  as  denning  numerical  convergence 
for  all  problems. 

Results  were  obtained  for  the  basic  iteration;  all  four  variations  of  the  VA 
method;  the  Ve  method;  orders  one,  two,  and  three  for  the  MVe  method;  and  for  k 
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values  of  two  through  five  (or  the  dimension  of  the  problem  if  less  than  five)  for  the 
RRE,  MPE,  and  AND  methods.  In  addition,  the  AND  method  was  applied  with 
the  relaxation  option.  However,  this  method  never  obtained  results  that  matched 
those  of  Anderson’s  higher  degree  methods;  hence,  results  of  the  relaxation  method 
will  be  shown  for  Example  1  only.  In  addition,  the  convergence  of  the  VA,  Vs, 
and  MVe  methods  were  usually  much  slower  than  that  of  the  AND,  MPE,  or  RRE 
methods.  Therefore,  the  results  of  the  VA,  Vs,  and  MVj  methods  will  not  be  shown 
for  examples  where  their  results  did  not  compare  favorably  with  the  other  methods. 
Because  of  the  structure  of  the  Vs  method,  the  results  will  show  the  difference 
vector  of  £°_15  if  n  is  even,  or  sjl_1,  if  n  is  odd.  The  RRE,  MPE,  and  AND  methods 
consistently  obtained  the  best  results,  faster  convergence.  For  each  of  these  three 
methods,  the  results  for  the  value  of  k  which  obtained  fastest  convergence  for  that 
particular  method  will  be  shown.  In  addition,  some  examples  will  show  results  for 
varying  values  of  k  for  a  particular  method. 

Before  we  continue,  it  should  be  understood  that  throughout  this  study  we 
have  assumed  that  the  fewest  number  of  iterations  implies  the  best  results.  This 
is  not  always  the  case.  Because  of  the  sophistication  of  some  of  the  methods  and 
the  simplicity  of  the  basic  iteration  scheme  of  some  problems,  more  extrapolations 
and  fewer  iterations  may  be  more  time  consuming  than  the  convergence  of  the  basic 
iteration.  However,  this  type  of  problem  is  not  in  the  majority,  and  fewer  iterations 
usually  means  less  computer  time  and  better  results. 

EXAMPLE  1:  The  first  example  is  a  simple  highly  degenerate,  linear  problem 
(Anderson,  1965,  Eq.  (5.1)).  Define  F(x)  of  Equation  (1)  by 

F(x )  =  Ax  —  <£6  =  0, 
where  d  is  a  free  parameter  and 
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Oij  =  < 


d, 

1, 


if  i  =  j 
if  i  ±  j. 


Choose  b  such  that  the  component  elements  of  the  solution  s  =  (s1} . . . ,  im)r  is 
Si  —  2/z,  i  =  1 ,m.  The  matrix  A  can  be  written  blsA  =  L  +  UP  +  D ,  where 
L  and  UP  are  defined  by  (86)  and  (87),  respectively;  and  D  is  the  diagonal  matrix 
[on, . . .  ,amm].  Therefore,  the  problem  can  be  rewritten  in  the  Jacobian  iteration 


form 


x  =  G(x)  =  -D~1{L  +  UP)x  +  D~xdb. 


(91) 


For  m  =  20, d  =  25,  and  xq  =  (l,...,l)r;  the  basic  iteration  sequence  con¬ 
verges  in  fifteen  iterations.  All  six  acceleration  methods  also  obtain  convergence 
(results  are  not  shown).  The  RRE,  MPE,  and  AND  methods  converge  after  only 
four  iterations  ( k  =  2).  The  Ve  and  MVe  methods  converge  in  5  iterations,  while 
the  VA  method  requires  18  iterations. 

When  d  is  set  to  15,  the  problem  is  a  little  different  because  the  basic  iteration 
scheme  produces  a  divergent  sequence.  However,  all  acceleration  methods  obtain 
the  correct  solution.  Figure  20  (page  95)  shows  best  results  for  all  acceleration 
methods  except  the  VA  method.  Figure  21  (page  95)  shows  results  of  the  basic 
iteration,  the  VA  method,  the  relaxation  option  of  the  AND  method  (relax),  and 
the  MPE  method  for  k  =  3  and  4.  Best  results  axe  for  order  2  for  the  MVe  method 
and  k  =  2  for  the  three  methods  using  k  values.  As  mentioned  earlier,  the  VA 
method  is  inappropriate  for  divergent  problems  unless  the  vectors  xn  and  xn+2  are 
interchanged  in  the  formula.  The  graph  of  the  VA  method  as  shown  in  Figure  21  is 
without  this  change  and  is  shown  for  comparison  purposes  only.  The  convergence 
rate  of  the  other  five  methods  are  similar. 

One  point  of  interest  noted  in  Figure  20  is  that  the  AND  method  converges 
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Figure  20.  Results  for  Example  1:  AND,  MPE, 
RRE,  MVe,  and  Ve 


ITERATION 

Figure  21.  Results  for  Example  1:  VA,  MPE 
(k  =  3  and  4),  and  Relax 
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for  k  >  2  in  four  iterations  due  to  its  ability  to  extrapolate  after  only  two  iterations 
even  though  the  final  value  of  k  may  be  larger  than  2.  For  the  RRE  and  MPE 
methods,  the  number  of  iterations  increases  for  k  >  2  is  shown  in  Figure  21  for  the 
MPE  method  with  k  =  3  and  4. 

EXAMPLE  2:  Consider  the  divergent  linear  problem  (91)  (Smith,  Ford,  and 
Sidi,  1987).  The  iteration  scheme  for  this  problem  was  discussed  in  Chapter  X, 
and  the  results  are  shown  in  Figure  22  for  the  RRE,  MPE,  AND,  and  MVe  meth¬ 
ods.  Even  though  the  dimension  of  this  problem  is  four,  the  degree  of  the  minimal 
polynomial  is  three  since  the  matrix  A  has  one  zero  eigenvalue.  All  three  methods 
are  able  to  detect  the  zero  eigenvalue  as  demonstrated  by  the  fact  that  the  best 
results  are  obtained  for  k  =  3  when  the  first  iteration  is  discarded.  In  addition,  the 
AND  method  achieves  machine  accuracy  (the  norm  of  the  difference  vector  equals 


Figure  22.  Results  for  Example  2:  AND, 
MPE,  RRE,  and  MVe 


zero)  for  k  =  4  and  the  first  iterate  discarded.  However,  the  number  of  iterations 
remains  seven  and  the  accuracy  of  the  sixth  iterate  is  not  as  high  as  the  sixth  iterate 
for  k  =  3.  If  the  first  iteration  is  not  discarded,  the  AND  method  still  converges  in 
seven  iterations  while  the  RRE  and  MPE  methods  require  one  additional  iteration 
each  to  converge.  However,  for  all  three  methods,  the  best  results  are  obtained 
with  k  =  4.  The  results  of  the  MVs  method  are  for  order  3  with  the  first  iterate 
discarded. 

There  is  another  point  of  interest  concerning  this  example.  Since  the  problem 
is  linear  and  the  degree  of  the  minimal  polynomial  is  three  ( k  =  3),  the  exact  solution 
should  be  found  by  all  three  methods  in  k  +  1  =  4  iterations  and  one  extrapolation. 
This  was  not  done  because  computer  computation  is  not  exact  arithmetic;  hence, 
rounding  errors  result  and  exact  convergence  in  four  iterations  is  not  obtained. 

Smith,  Ford,  and  Sidi  (1987)  obtained  a  “wrong”  solution  for  this  problem 
using  the  RRE  method.  The  reason  their  RRE  method  failed  to  converge  to  the 
correct  solution  is  discussed  in  the  Appendix. 

EXAMPLE  3:  For  a  final  look  at  a  linear  case,  we  find  the  solution  of  the 
system 

5  7  6  5  23 

7  10  8  7  32 

x  = 

6  8  10  9  33 

5  7  9  10  31 

J 

by  use  of  the  Jacobian  iteration,  (91),  with  an  initial  vector  x0  =  (0, 0, 0, 0)r  (Smith, 
Ford,  Sidi,  1987,  Ex.  1).  The  solution  is  (1,1, 1,1).  One  of  the  eigenvalues  for  this 
problem  is  near  0.9985,  causing  the  matrix  I  -  A  to  be  nearly  singular.  Figure 
23  (page  98)  gives  the  results  for  the  RRE  and  MPE  methods  (k  =  4),  the  AND 
method  (k  =  4  and  3),  and  the  MVe  (order  3).  This  is  a  linear  problem  whose 
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Figure  23.  Results  for  Example  3:  AND, 
MPE,  RRE,  and  MVe 


minimal  polynomial  has  degree  four.  Therefore,  ignoring  roundoff,  convergence 
should  be  obtained  in  k  +  1  =  5  iterations  and  one  extrapolation.  However,  once 
again,  this  is  not  achieved  due  to  rounding  errors. 

One  observation  is  that  whereas  the  AND  methods  converge  in  7  and  11 
iterations,  the  RRE  and  MPE  methods  do  not  converge  to  high  accuracy  in  even  25 
iterations.  In  fact,  the  convergence  of  these  two  methods  does  not  improve  in  250 
iterations  for  any  value  of  k.  For  both  methods,  the  log  of  the  difference  vectors 
fluctuates  between  —9  and  —13  with  no  set  pattern.  As  a  result,  the  log  of  the 
error  vectors  never  gets  smaller  than  —11.  This  inability  to  improve  convergence 
beyond  a  certain  value  is  referred  to  as  “limited  accuracy.”  The  VA  method  and  the 
e  methods  have  the  same  difficulty  on  this  problem  only  with  much  lower  accuracy. 

There  were  several  variations  made  to  the  basic  iteration  scheme  to  attempt 
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to  achieve  convergence.  These  included  using  combinations  of  different  relaxation 
factors,  discarding  the  first  few  iterates,  and  discarding  a  set  number  of  iterates  at 
the  beginning  of  each  extrapolation  throughout  the  run  of  the  program  or  only  after 
the  limited  accuracy  had  been  achieved.  However,  no  improvement  was  made  as  all 
attempts  eventually  led  to  some  value  of  limited  accuracy. 

Before  investigating  the  cause  of  this  limited  accuracy  problem,  a  few  defi¬ 
nitions  are  in  order.  Given  the  vector  norm  ||xj|,  define  the  norm,  ||A|j,  and  the 
condition  number,  cond(A),  of  the  matrix  A  by 

Mil  =  S»P  w  cond(4)  =  ||A||||A-‘||, 

**d  If  II 

— ^ 

respectively  (Conte  and  de  Boor,  1980).  For  the  linear  system  Ax  =  6,  define  the 
relative  error  as  where  x  is  the  computed  vector  of  x,  and  the  residual  error 

by  It  is  shown  by  Golub  and  Van  Loan  (1983)  how  the  condition  number 

- — ♦ 

and  the  relative  errors  in  A ,  x,  and  6  relate.  If  cond(A)  «  1,  then  the  relative  error 
and  the  residual  error  will  be  of  the  same  order  of  magnitude;  hence,  x  is  a  good 
approximation  of  x.  If  the  condition  number  is  large,  then  a  small  change  in  the 
data  MAY  cause  a  large  change  in  the  solution.  In  short,  the  condition  number 
“quantifies  the  sensitivity  of  the  Ax  =  b  problem”  (Golub  and  Van  Loan,  1983). 
Ortega  and  Pt'  le  (1981)  give  the  example  that  a  condition  number  of  10®  could 
result  in  a  loss  of  6  decimal  digits  of  accuracy.  A  matrix  with  a  large  condition 
number  is  called  an  “ill-conditioned”  matrix. 

The  RRE,  MPE,  and  AND  methods  all  involve  solving  a  linear  system  Ux  =  b 
to  determine  the  extrapolated  vector.  Therefore,  consider  the  condition  number 
of  U  for  each  extrapolation.  For  the  k  values  shown  in  the  Figure  22,  the  first 
condition  number  for  the  RRE  method  is  1013  and  for  the  MPE  method  is  lO10. 
The  remaining  condition  numbers  vary  between  103  and  109.  For  the  AND  method, 
the  first  condition  number  is  104  and  the  condition  number  for  the  matrix  when 
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convergence  is  obtained  is  2.67.  This  is  a  major  difference  and  could  be  a  factor  in 
allowing  Anderson’s  method  to  overcome  the  limited  accuracy  problem. 

Since  the  example  is  a  small  linear  problem,  it  can  be  solved  by  Gaussian  elim¬ 
ination.  The  condition  number  of  the  matrix  A  is  50  and  the  computed  solution 
has  a  relative  error  of  0.1896  x  10-13.  The  relative  error  of  Anderson’s  computed 
solution  is  0.2251  x  10~12.  The  relative  errors  for  the  RRE  and  MPE  methods  never 
obtain  an  accuracy  higher  than  0.4467  x  10~10  and  0.1973  x  10~8,  respectively.  The 
AND  method,  even  with  the  low  condition  numbers,  still  cannot  achieve  quite  the 
same  degree  of  accuracy  that  the  Gaussian  elimination  method  achieves.  However, 
when  comparing  acceleration  methods,  the  AND  method  definitely  shows  superi¬ 
ority  in  overcoming  the  problem  of  limited  accuracy.  The  RRE  method  for  Smith, 
Ford,  and  Sidi  (1987)  failed  to  converge  for  this  problem.  See  the  Appendix  for 
more  details. 

EXAMPLE  4:  The  first  nonlinear  example  is  a  quadratic  problem  with  solu¬ 
tions,  (1,1, 1,1)  and  (3, 3, 3, 3)  (Gekeler,  1972,  Ex.  V).  Define  G(x)  of  equation  (2) 
by 

G(x)  =  Ax  +  b  +  Q{x), 

where 


A  = 


3.9 

-3.7 

2.4 

-0.6 

2.4 

-2.0 

2.2 

-0.6 

2.4 

-3.6 

4.1 

-0.9 

J 

2.8 

-5.2 

4.8 

-0.4 

-0.75(1,1,1, 

If, 

Q(x) 

=  — 0.25( 

b  : 

£0  =  1.5(1, 1,1, l)r. 


.2  \T 


The  basic  iteration  converges  to  the  solution  (3, 3, 3, 3)  in  52  iterations  with  a  con- 
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vergence  criterion  of  1  x  10-14,  C14;  however,  in  400  iterations  the  norm  of  the 
difference  vector  never  gets  smaller  than  the  normal  convergence  criterion,  C15. 

This  problem  has  some  unusual  properties.  Figure  24  (page  102)  shows  the 
results  of  the  basic  iteration,  the  Ve  method,  the  MVe  method  (order  4),  and  the 
MPE  and  RRE  methods  with  k  =  4.  Each  of  these  methods  converges  to  (3, 3, 3, 3). 
Figure  25  (page  102)  shows  the  convergence  of  the  VA  with  Q  equal  to  formula  (33), 
the  AND  method  ( k  =  3),  the  MVe  (order  3),  and  the  MPE  method  (k  =  1).  These 
four  methods  converge  to  the  solution  (1,1, 1,1).  In  addition,  Figure  25  shows  the 
convergence  of  Newton’s  iterative  method,  the  “most  famous  iterative  method  for 
obtaining  roots  of  equations  (as  well  as  for  solving  systems  of  nonlinear  equations 
...).”  (Ortega  and  Poole,  1981,  p.  128).  Convergence  is  not  obtained  for  the  RRE 
method  ( k  <  4)  and  the  MPE  method  (k  =  2  and  3).  For  these  methods,  the 
system  either  overflows  or  the  limited  accuracy  is  10"1. 

A  possible  reason  for  this  difference  for  the  methods  involving  k  values  is  the 
matrix  determined  for  computing  the  first  extrapolated  vector.  For  both  methods 
that  converge  to  (3,  3, 3, 3),  the  matrix  is  singular.  As  a  result,  the  matrix  of  the 
linear  system  has  a  rank  smaller  than  k ;  hence,  a  linear  system  of  smaller  degree  is 
solved  resulting  in  a  first  computed  extrapolated  vector  near  (3,  3, 3, 3).  Remaining 
iterates  and  extrapolations  converge  to  this  vector.  The  rank  of  the  matrix  for 
the  other  converging  methods  is  k  and  the  resulting  extrapolated  vector  is  near 
(1,1, 1,1). 

As  stated  in  Chapter  I,  small  nonlinear  problems  will  in  practice  be  solved 
by  methods  other  than  those  discussed  in  this  thesis.  To  illustrate  this  point,  this 
problem  is  solved  by  Newton’s  method.  Let  F(x )  =  G(x )  -  x,  then  Newton’s 
method  is  the  iterative  formula 


xn+1  =  xn  +  y„,  where  F'(xn)yn  =  -F(xn) 
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Figure  24.  Results  for  Example  4:  Methods 
Converging  to  (3,3,  3,  3) 
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for  F'(xn)  the  Jacobian  matrix  of  F(x)  evaluated  at  xn.  As  can  be  seen  in  Figure 
25,  all  the  methods  that  converge  to  (1,1, 1,1)  have  some  difficulty  with  limited 
accuracy.  However,  Newton’s  method  and  the  AND  method  overcome  this  problem 
and  achieve  convergence  fairly  quickly. 

Example  5:  The  next  example  comes  from  Wynn  (1964,  Eq.  (1))  and  is 
referred  to  as  the  Lichtenstein-Gershgorin  integral  equation.  The  iteration  scheme 
for  n  =  0, 1, . . .  is 

e  fxi  =  *r  _*&(!> _ Mn(7r-Q  , 

n+1  IT  JO  1  —  JC2  COs(t  +  x)  1  +  &2  COs(<  +  x) 

_  k~x  sin(x) 

+2  arctan  r— - - — - - . 

[1  —  cos(x)]  [fc3  cos(x)  —  A:-2] 

where 


k\  —  ( k 2  +  1)  1,  &2  =  k\(k2  —  1),  fc3  =  1  —  k  2,  and  9q  =  0. 

The  integrals  are  approximated  by  the  trapezoid  rule  with  end  corrections: 

/a+mh  f  l  1  -i 

f{t)  dt  =  k  |  -f0  +  f\  H - V  fm-1  +  ~fm  +  Cj  (92) 

for 

1  1  IQ  1 

C  =  ^(A/o  -  v/n)  -  — (A2/0  +  V2/n)  +  ^(A3/o  -  V3/n)  -  ^(AVo  +  VVn), 

where 

A/0  =  h  ~  fo,  An/0  =  An+1/0  -  A  Vo, 

V/n  =  fn—l  -  /»,  and  Wn  =  V""7n  -  VVn. 

Choosing  m  =  73  and  x<j  =  (0, . . . ,  0)r,  the  basic  iteration  converges  in  136 
iterations.  Figure  26  (page  104)  shows  the  graphs  for  the  convergence  of  the  RRE, 
MPE,  and  AND  methods  with  k  =  5,  the  MVf  method  (order  3),  and  the  VA 
method  with  Q  equal  to  formula  (33).  As  can  be  seen  from  the  figure,  this  example 
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Figure  26.  Results  for  Example  5:  AND, 
MPE,  RRE,  VA,  and  MVe 


has  very  similar  results.  This  problem  does  illustrate  the  importance  of  an  accel¬ 
eration  method.  Due  to  the  nature  of  an  integral  equation  problem,  the  computer 
time  required  to  compute  one  iterate  is  much  more  than  any  of  the  previous  ex¬ 
amples.  Therefore,  to  decrease  the  number  of  iterations  to  less  than  one-fifth  that 
required  by  the  basic  iteration  is  a  significant  reduction  in  computer  time  and  cost. 
There  are  no  other  interesting  points  concerning  the  results  of  this  example  other 
than  the  fact  that  this  is  the  first  example  in  which  a  method  other  than  the  AND 
method  even  came  close  to  having  the  best  results.  The  AND  method  is  still  quite 
competitive,  however. 

Example  6:  The  next  examples  are  the  integral  equation  problems  (Anderson, 
1965,  Eqs.  (5.10)  k  (5.11)): 

/*(*)  =  f  x  f{t)  (cos  *  ■)  dt  -  J  and  (93) 


J 
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„  .  3?r\/2  y1  ,,,  .  7r|x  —  t\  ,  1  7TX 

f(x )  =  — ^r-  y_1  f (0 cos — 4 — ^~4cosX’  (94) 

with  both  solutions  /(x)  =  cos(ttx/4). 

Letting  m  =  101, /o  =  (l,...,l)r,  and  integrating  both  problems  with  the 
iteration  scheme  (92);  Integral  (93)  converges  in  75  iterations  while  Integral  (94) 
diverges.  Results  for  Integral  (93)  axe  shown  in  Figure  27  (page  106).  The  AND 
method  clearly  obtains  the  fastest  convergence.  Once  again,  the  RRE  and  MPE 
methods  ( k  =  2)  have  similar  convergence.  This  example  will  also  show  how  using 
a  convergent  generated  sequence  in  reverse  order  can  affect  the  convergence  rate.  If 
the  MPE  method  ( k  =  2)  is  applied  to  the  sequence  S  =  {xo, . . . ,  Xfc+1},  the  norm 
of  the  thirteenth  difference  vector  is  10~u  as  compared  to  10-15  when  the  method 
is  applied  to  the  sequence  T  =  {x*+1,...,x o}.  In  addition,  it  requires  four  more 
iterations  to  converge  to  10-15.  Though  there  are  exceptions  to  the  riile,  applying 
the  acceleration  method  to  the  sequence  T  instead  of  the  sequence  S  usually  reduces 
the  number  of  iterations  for  convergence  if  5  is  a  convergent  sequence.  Using  the 
most  accurate  estimate  of  the  solution  as  the  first  term  of  the  sequence  will  usually 
cause  faster  convergence. 

For  Integral  (94),  the  basic  iteration  scheme  produces  a  sequence  that  diverges 
quadratically,  the  logarithm  of  the  norm  of  the  difference  vector  roughly  doubling 
as  n  increases  by  one.  Unlike  many  linear  divergences,  this  divergence  cannot  be 
“tamed”  by  using  a  relaxation  factor  in  the  interval  (0, 1).  The  two  e  methods  also 
produce  divergent  sequences.  Because  the  norm  of  the  difference  vectors  increases 
rapidly  even  for  the  first  few  iterates,  the  RRE  and  MPE  methods  for  k  >  2  do 
not  obtain  convergence.  For  k  =  1,  both  methods  converge,  but  the  convergence  is 
slow  as  illustrated  in  Figure  28  (page  106).  The  graphs  clearly  show  that  the  AND 
method  gives  the  best  results.  In  fact,  for  all  values  of  the  results  of  the  AND 
method  are  almost  identical. 
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Example  7:  This  problem  is  the  “model  problem”  of  Varga  (1962)  and  Young 
(1971).  It  is  designed  to  solve  the  Dirichlet  problem  on  a  rectangle.  Define  fl  as  the 
interior  of  a  rectangle,  S  as  the  boundary  of  the  rectangle,  and  flu  S  as  their  union. 
Let  G(x,y)  and  y(x,y)  be  continuous  functions  defined  on  R  and  S,  respectively. 
Then  the  desired  solution  is  a  function  u (x,y)  that  is  continuous  on  flU  S,  is  twice 
continuously  differentiable  on  fl,  and  satisfies  Poisson’s  equation 


d2u  d2u 

a*  +  w  =  G{x’y)'  (95) 

In  addition,  u(x,y)  =  g(x,y)  on  S.  If  G(x,y )  =  0,  then  (95)  reduces  to  Laplace’s 
equation. 

The  function  u  is  found  numerically  by  finding  approximations  to  the  function 


at  a  finite  number  of  interior  points.  These  points  are  obtained  by  superimposing 
a  rectangular  mesh  of  horizontal  and  vertical  lines  with  uniform  spacing.  With 
reference  to  Figure  29  (page  108),  define  (x0,yo)  and  ( xp,ym ),  p  and  m  integers, 
as  the  lower  left  point  and  upper  right  point,  respectively,  of  the  rectangle.  Also 
define  d  —  xp  —  xo,w  =  ym  —  yo,h  =  d/p,  and  k  =  w/m.  Then  the  spacing  of 
the  rectangular  mesh  is  h  for  the  vertical  lines  and  k  for  the  horizontal  lines.  The 


spacing  in  Figure  29  is  h  =  d/Z  and  k  =  w/4.  Other  points  of  the  mesh  are 
(xiiVj)  =  (*o  +  hi, y0  +  kj).  Denote  the  functional  value  u(xi,y_7),z  =  0, . . .  ,n  and 
j  =  0,. . .  ,m,  by  Uij  =  u(xj,t/J).  Hence  the  solution  will  be  the  approximations  U{j. 

Finite  difference  approximations  to  the  second  derivative  with  respect  to  x 
and  y  are  defined  by 

0  =  [u(x  +  h,y)~  2u(x,  y)  +  u(x  -  h,  y)]  /h2  and 

(96) 

0  =  [«(*»y  +  *)  ~  2 i*(x,y)  +  u(x,y-  *)]  /k2, 
respectively  (Ortega  and  Poole,  1981).  To  simplify  the  problem,  assume  the  re¬ 
gion  i*U  5  is  the  unit  square,  h  =  fc,(x0,y0)  =  (0,0),  and  G(x,y)  =  y(x,y)  =  0. 
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Figure  29.  Rectangular  Mesh  with  Spacing  of 
h  =  d/3  and  k  =  w/4 


Therefore,  adding  the  two  equations  of  (96)  and  using  (95)  result  in  a  natural 
iteration  formula  for  the  problem: 


u 


("+i)  _ 


•j 


-  L(n)  .  +  «(")  +  „(n)  +  u(")  1  /4 

—  +  U»-l  J  +  ui,j+ 1  + 


(97) 


i,j  =  1.  The  right-hand  side  of  (97)  is  referred  to  as  the  five-point  Jacobian 

operator. 

The  basic  iteration  generates  a  slowly  convergent  sequence  as  illustrated  by 
the  fact  that  the  infinity  norm  of  the  272nd  difference  vector  is  only  0.944519  x  10-7 
to  six  place  accuracy.  Figure  30  (page  109)  shows  the  first  50  iterations  of  the  best 
results  of  the  RRE,  MPE,  and  AND  methods.  The  other  acceleration  methods  do 
not  work  well  on  this  problem.  As  an  example,  all  three  orders  of  the  MVe  method 
converge  slower  than  or  similar  to  the  basic  iteration,  as  shown  in  the  figure  for 
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Figure  30.  Results  for  Example  7:  AND, 
MPE,  RRE,  and  MVe 


order  three.  This  problem  is  a  good  example  to  illustrate  how  using  a  relaxation 
factor  other  than  unity  and  discarding  a  few  iterations  may  provide  a  faster  conver¬ 
gence.  Testing  was  done  for  relaxation  factors  between  0.5  and  1.5  and  discarding 
one  to  five  iterations.  Though  not  all  variations  produced  faster  convergence,  better 
results  than  those  shown  in  Figure  30  were  obtained  for  all  three  methods.  Table 
15  (page  110)  gives  some  of  the  results  obtained  and  the  variations  for  the  AND 
and  RRE  methods.  It  should  be  noted  that  this  is  only  the  second  test  problem  in 
which  another  method  matched  the  results  of  the  AND  method. 

Example  8:  This  example  is  another  nonlinear  integral  equation  (Rail,  1969) 

F(x)  =  1  -I-  (tt0/2)xF(x)  f  F(y)_  jy,  0  <  x  <  1,  (98) 

Jo  x  -f  y 

for  0  <  7r0  <  1.  The  background  material  for  the  equation  is  fairly  elaborate  (Chan¬ 
drasekhar,  1960).  Due  to  the  natural  iterative  form  of  (98),  the  basic  iteration 


no 


TABLE  15 

ITERATIONS  REQUIRED  FOR  CONVERGENCE  OF 
EXAMPLE  7  WITH  MODIFICATIONS  TO  THE 
ANDERSON  AND  RRE  METHODS  WITH  A 
CONVERGENCE  CRITERION  OF  10"7 


Method 

Relaxation 

Parameter 

Number  of  Discarded 
Iterations 

Number  of 
Iterations 

AND 

0.5 

2 

37 

AND 

1.0 

1 

45 

RRE 

0.5 

1 

45 

RRE 

1.0 

1 

44 

formula  is 

Fn+ i(*)  =  1  +  (iro/2 KW  /'  dy,  (99) 

Jo  x  +  y 

n  =  0,1,...,  with  Fo(x)  =  1,  0  <  x  <  1.  However,  analytic  difficulties  do  develop 
involving  the  integration  portion  for  finding  F2(x).  Rail  showed  that  for  F0(x)  =  1, 
for  all  x  an  element  of  the  interval  of  integration,  to  be  a  satisfactory  initial  approx¬ 
imation  to  the  solution,  tw0  is  restricted  by  0  <  tt0  <  (-v/2  —  l)/(ln2)  ~  0.59758 . . . 
Therefore,  he  constructed  a  corresponding  arithmetic  model  by  introducing  a  “nu¬ 
merical  integration  ride”  of  the  form 

JQ  f{s)  ds  ->  f^Wi  f{sx ), 

where  s;,  0  <  s,  <1,  i  =  1, . . .  ,m,  axe  nodes;  the  parameters  wx,  i  =  1, . . .  ,m,  are 
weights;  and  m  is  the  order  of  the  rule.  The  integral  portion  of  (98)  becomes 


(100) 


Ill 


The  solution  F(x )  is  approximated  by  determining  F(xi ),  0  <  x*  <  1, 
t  =  1, . . .  ,m.  Choosing  Xi  =  y*,  t  =  1, . . .  ,m  and  using  (100),  the  value  at  x<  is 

T»  .  in  . 

*•(*0  =  1  +  (iro/2)F(*0  £  Fix/). 

j=l 

Letting  6^  =  *'+J. ,  Equation  (98)  becomes 

m 

/i  =  1  +  (tto/2)/,-  i  =  1, . . . , m, 

»=i 

where  /;  =  F(xj).  Defining  the  m-dimensional  vectors  z  and  1  by  x  =  (/i,...,/m)r 
and  1  =  (1,...,1)T,  respectively,  the  iteration  formula  takes  the  form 

zn+x  =  14-  (7r0/2)zn  (g)  Fxn, 

where  B  =  [&ij]  and  ®  stands  for  the  component-by-component  multiplication  of 
the  vectors:  y®z  =  {yxzl, . . .  ,ymzm)  for  y  =  (yi, . . .  ,ym)r  and  z  =  (z1? . . . ,  zm)r. 

Table  16  gives  the  nodes  and  weights  to  seven  places  for  the  Gaussian  integra¬ 
tion  rule  of  order  nine  (Milne,  1949).  Using  these  values  and  7r0  =  0.1(i),  for 
i  =  1 , . . . ,  10,  Rail  obtained  convergence  to  eight  decimal  places  for  all  cases; 


TABLE  16 

NODES  AND  WEIGHTS  FOR  THE  GAUSSIAN 
INTEGRATION  RULE  OF  ORDER  NINE 


i 

Si 

Wi 

t 

Si 

Wi 

1 

0.0159199 

0.0406372 

6 

0.6621267 

0.1561735 

2 

0.0819844 

0.0903241 

7 

0.8066857 

0.1303053 

3 

0.1933143 

0.1303053 

8 

0.9180156 

0.0903241 

4 

0.3378733 

0.1561735 

9 

0.9840801 

0.0406372 

5 

0.5000000 

0.1651197 
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TABLE  17 

NUMBER  OF  ITERATIONS  REQUIRED  TO  OBTAIN 
CONVERGENCE  ON  RALL’S  PROBLEM  FOR 
DIFFERENT  VALUES  OF  tt0 


7T  o 

Number  of 
Iterations 

*0 

Number  of 
Iterations 

0.1 

7 

0.6 

17 

0.2 

8 

0.7 

21 

0.3 

10 

0.8 

28 

0.4 

12 

0.9 

43 

0.5 

14 

1.0 

10587 

however,  as  7r0  increased  so  did  the  number  of  iterations  for  convergence,  Table 
17  (Rail,  1969).  Table  17  also  shows  that  for  ir0  =  1,  convergence  is  extremely 
slow,  and  according  to  Rail  has  limited  accuracy.  Because  w0  =  1  is  by  far  the  most 
difficult  case,  the  acceleration  methods  are  applied  to  this  problem  for  this  case  only. 
Figure  31  (page  113)  shows  results  obtained  for  the  MPE,  RRE,  and  AND  methods 
for  the  first  50  iterations.  What  the  figure  does  not  show  is  that  this  problem 
is  another  example  of  limited  accuracy  for  the  RRE  and  MPE  methods.  Both 
methods  converge  to  CIO,  but  neither  one  converges  to  Cll  in  3000  iterations  even 
with  modifications  to  the  relaxation  factor  and  the  number  of  iterations  discarded. 
It  should  be  added  that  the  AND  method  converges  to  C15  for  all  values  of  k  >  2 
in  less  them  132  iterations. 

Example  9:  The  next  example  has  the  largest  dimensional  value  of  all  the  test 
problems.  It  approximates  the  steady-state  solution  of  the  scalar  three-dimensional 
Burger’s  equation 


ut  +  u(ux  +  uv  +  ut)  =  eAu, 


(101) 
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Figure  31.  Results  for  Example  8:  Rail’s 
Problem  with  7r0  =  1 


on  the  unit  cube  (Hyman  and  Manteuffel,  1984).  The  parameter  t  represents  time; 
ut,ux,uy,  and  uz  are  the  partials  of  u  with  respect  to  f,x,y,  and  z,  respectively; 
and  A  is  the  Laplacian  operator 

d2u  d2u  d2u 

d^+d^+d^’ 

Hyman  and  Manteuffel  tested  an  acceleration  method  they  developed  by  studying 
the  convergence  of  a  second-order  Runge-Kutta  method  for  this  problem  using 
e  =  0.02  and  the  Dirichlet  boundary  conditions 

u(0,  y,  z)  =  u(x,  0,  z)  =  ti(x,  y,  0)  —  0  and 

=  **(*,!,*)  =  «(*,y,i)  =  1 

to  provide  a  thin  boundary  layer.  In  addition,  a  time  step  of  A t  =  0.5(Ax)  was 
chosen.  As  in  Example  7,  a  second  order  finite  difference  equation  was  used  to 
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approximate  the  solution  on  a  uniform  grid  of  N  points.  They  gave  results  for 
N  =  8000  (a  20  x  20  x  20  grid  of  points)  (Hyman  and  Manteuffel,  1984,  p.  312). 

The  exact  structure  of  their  test  case  is  not  cleax.  First,  there  is  a  class  of 
second  order  Runge-Kutta  integration  methods,  and  Hyman  and  Manteuffel  do  not 
mention  which  one  of  these  methods  they  used.  Gear  (1971)  defined  the  two-step 
calculation  of  a  Runge-Kutta  method  as 

9i  =  yn  +  ahf(yn,tn) 

Vn+ 1  =  yn  +  bhf(yn,tn)  +  chf(qlytn  +  dh), 

where  dy/dt  =  f(y,t)’,  a,b,c,  and  d  are  parameters;  and  h  =  Ax.  To  make  the 
expansion  of  yn+i  and  the  Taylor  series  agree  as  closely  as  possible,  the  relationship 
between  the  parameters  must  be  b  =  1  —  c  and  a  =  d  =  c/2.  Therefore,  the  basic 
iteration  formula  for  Burger’s  equation  with  u  =  y,  A t  =  0.5/i,  and  ut  =  /(y,  t )  is 

9i  =  un  +  a(At)(ut)n 
tin+i  =  un  +  6(At)(ut)  +  c(Ai)(ut)n+1, 

where  (ut)n+i  is  ut  evaluated  at  qi  and  tn  +  dh.  Three  common  second  order  Runge- 
Kutta  methods  are  for  c  =  1/2,  3/4,  and  1  (Gear,  1971,  p.  31).  A  second  unclear 
area  is  whether  the  number  of  grid  points,  N,  includes  the  boundary  points.  Because 
the  software  and  the  exact  parameters  for  their  test  problem  were  not  available,  the 
results  shown  for  this  example  are  for  N  =  8000  to  be  the  number  of  interior  points 
and  for  c  =  1/2.  Results  of  the  basic  iteration  do  not  exactly  match  those  shown 
by  Hyman  and  Manteuffel;  however,  the  problem  is  still  a  good  test  problem  due 
to  the  size  of  its  dimension. 

Results  obtained  for  this  problem  are  shown  in  Figure  32  (page  115)  for  the 
AND  ( k  —  4),  MPE  (fc  =  3),  RRE  (k  =  1),  and  MVe  (order  3)  methods  for  the 
first  50  iterations.  In  addition,  the  graphs  plot  the  infinity  norm  of  the  error  vector, 
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Figure  32.  Results  for  Example  9:  Hyman  and 
Manteuffel’s  Test  Problem 


e*n-i  =  i*— x„_i,  instead  of  the  difference  vector.  The  change  is  to  match  Hyman  and 
Manteuffel’s  test  problem  as  close  as  possible.  The  only  methods  that  consistently 
increase  the  convergence  rate  are  the  AND  and  MPE  methods.  The  other  methods 
do  not  work  very  well  for  this  problem.  In  fact,  the  convergence  of  the  MPE  and 
RRE  methods  is  not  smooth  as  illustrated  by  the  periodic  peaks  in  their  graphs 
even  though  the  basic  iteration  generates  a  convergent  sequence.  In  addition,  the 
MPE  method  is  very  erratic.  Hence,  these  two  extrapolation  methods  have  their 
problems  in  this  example. 

Example  10:  The  last  example  comes  from  Moler  (1967).  The  problem  is 
to  solve  for  the  eigenvalues  and  eigenfunctions  of  the  Laplacian  operator,  Equa¬ 
tion  (102)  in  two  variables  only,  on  an  “L”  shaped  region  L ,  Figure  33  (page  116). 


116 


The  eigenvalues  q  and  the  functions  u(p),  not  identically  zero,  are  to  satisfy 


A u(p)  +  qu(p)  =  0,  p  =  (z,y)  an  element  of  L 
u (p)  =  0,p  an  element  of  L. 


(103) 


Since  there  are  infinitely  many  eigenvalues,  only  the  smallest  one  is  considered. 

Once  again,  the  solution  is  approximated  using  finite  differences  over  a  square 
mesh  of  width  h  =  1  /N,  N  an  integer.  Letting  =  u(xi>yj),  the  five-point 
Laplacian  operator  is  define  by 


^L^ij  —  "I"  Ui—lj  d"  d"  4Ujj]  /h  . 


(104) 


However,  direct  iteration  of  the  Laplacian  operator  will  produce  the  largest  eigen¬ 
value.  The  smallest  eigenvalue  can  be  found  by  use  of  the  five-point  Jacobian 
operator: 


AjUij  =  [ui+u  +  Ui-ij  +  Uij+i  +  /4. 


(105) 
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Figure  34.  Results  for  Example  10:  AND, 

MPE,  RRE,  MVe,  and  Ve 

Denote  qi  and  qj  as  the  eigenvalues  of  the  Laplacian  and  Jacobian  operators,  re¬ 
spectively.  Using  equations  (103),  (104),  and  (105),  the  following  relationship  holds: 

9l  =  [ui+ij  +  Ui-ij  -|-  mj+i  +  Uij-i  -  Auij]  lh7Uij 
=  4  +  ui- ij  +  1  +  «ij-i)/4«u  -  1]  /h* 

=  4 (qj  -  1  )/h*.  (106) 

Hence,  the  approach  of  finding  q,  the  smallest  eigenvalue,  is  to  solve  the  problem 
by  using  (105)  and  then  convert  the  solution  from  an  eigenvalue  of  the  Jacobian 
operator  to  one  of  the  Laplacian  operator  by  the  Relation  (106). 

Figure  34  shows  the  results  for  the  first  50  iterations  for  the  AND,  RRE, 
MPE,  and  both  e  methods.  For  this  last  problem,  the  graphs  plot  the  log  of  the 
Euclidean  norm  of  the  difference  vector  instead  of  the  infinity  norm.  The  RRE  and 
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MPE  methods  do  not  work  well  on  this  problem.  In  fact,  for  values  of  k  not  shown 
in  Figure  34,  the  convergence  is  slower  than  that  obtained  by  the  basic  iteration. 
The  AND  and  Ve  methods  do  accelerate  the  convergence,  but  even  the  accelerated 
convergence  is  slow.  However,  the  solution  can  be  obtained  in  fewer  iterations  by 
applying  either  one  of  these  methods.  (Note  that  the  Ve  method  has  enormous 
storage  and  time  requirements  in  this  example.) 


CHAPTER  XII 


THE  GENERALIZED  MINIMUM  RESIDUAL  ALGORITHM 

Information  on  another  acceleration  method,  the  Generalized  Minimum  Resid¬ 
ual  (GMRES)  algorithm,  was  received  just  prior  to  the  completion  of  this  thesis. 
Because  of  the  time  factor  and  the  complexity  of  the  software  of  the  method,  test¬ 
ing  of  the  method  was  minimal.  The  GMRES  algorithm  was  developed  for  linear 
systems  by  Saad  and  Schultz  (1986).  The  method  has  had  further  development 
by  Kerkhoven  and  Saad  (1987),  Brown  and  Saad  (1987),  and  Burkhart  and  Young 
(1988).  There  now  exist  routines  for  both  linear  and  nonlinear  problems.  Discussion 
of  this  method  will  be  sketchy. 

The  nonlinear  GMRES  method  resembles  Newton’s  method  for  solving  a  sys¬ 
tem  of  nonlinear  equations.  However,  GMRES  reduces  the  effective  dimension  of 
the  solution  space.  The  method  involves  finding  an  orthonormal  basis  of  the  Krylov 

subspace  Kk  =  span{vi ,  Avi , . . . ,  Afc_1,Ui},  where  A  is  the  matrix  of  the  linear  sys- 

— ♦  — ♦ 

tern  Ax  =  b  and  Vi  is  the  normal  vector  Vi  =  f0/f0 ,  f0  =  b  —  Ax, q.  The  basis  is 

found  by  a  procedure  called  Arnoldi’s  algorithm. 

Software  for  GMRES  was  obtained  from  Burkhart  and  Young  (1987)  of  Boe¬ 
ing  Computer  Services.  The  test  driver  program  that  was  provided  for  the  nonlin¬ 
ear  GMRES  routine  solved  Laplace’s  equation  on  a  square  mesh,  using  an  SSOR 
(symmetric  successive  over- relaxation)  iteration  (Varga,  1962)  with  a  mesh  size  of 
h  =  1/22  and  optimizing  the  over- relaxation  factor  OMEGA  automatically  dur¬ 
ing  the  process.  Nonlinear  GMRES  solved  this  problem  to  full  double  precision 
accuracy  using  a  total  of  129  sweeps  of  SSOR  over  the  grid. 
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Anderson’s  method  (fc  =  5)  on  this  problem,  wi;n  OMEGA  set  to  1.75,  which 
is  near  the  optimum  value,  requires  only  about  60  iterations  to  achieve  full  double 
precision  accuracy.  However,  the  GMRES  software  is  organized  in  a  very  conserva¬ 
tive  way  in  order  to  avoid  divergence.  If  similar  search  strategies  were  used  in  the 
Anderson  routine,  the  number  of  iteration  would  be  increased  considerably.  (Such 
search  strategies  should  be  an  option  of  the  software,  used  only  when  divergence  is 
anticipated  or  detected.) 


CHAPTER  XIII 


SUMMARY  AND  CONCLUSIONS 

The  intent  of  this  thesis  was  to  demonstrate  the  importance  of  acceleration 
methods  and  to  compare  several  of  these  methods  both  theoretically  and  numeri¬ 
cally.  For  each  method,  the  theory  and  the  algorithm  were  derived  for  the  linear 
case.  However,  through  numerically  testing  the  algorithms  on  different  types  of 
problems,  it  was  shown  that  the  methods  can  be  applied  to  both  linear  and  nonlin¬ 
ear  problems. 

Clearly,  the  purpose  of  acceleration  methods  is  to  reduce  the  number  of  it¬ 
erations  required  to  solve  numerically  a  mathematical  problem  in  a  vector  space. 
All  methods  presented  in  this  thesis  demonstrate  the  capability  of  achieving  this 
purpose,  though  the  convergence  rate  may  vary  for  different  problems.  This  in  itself 
is  of  great  value  since  for  the  majority  of  practical  problems  reducing  the  number 
of  iterations  also  reduces  the  computer  time  and  cost.  In  addition,  acceleration 
methods  also  have  demonstrated  the  capability  of  accelerating  some  divergent  se¬ 
quences  to  the  solution  of  a  problem.  Therefore,  a  greater  number  of  problems 
may  be  solved  numerically  by  applying  an  acceleration  technique  to  the  generated 
sequence. 

There  are  three  categories  of  acceleration  models:  the  static  model,  the  semi¬ 
dynamic  model,  and  the  fully  dynamic  model.  If  an  extrapolation  accelerates  the 
convergence,  then  one  may  suspect  that  the  fully  dynamic  model  will  provide  the 
fastest  convergence,  since  extrapolation  is  accomplished  after  every  iterate  once  the 
first  extrapolation  is  done.  The  only  fully  dynamic  method  for  vectors  presented, 
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Anderson’s  Generalized  Secant  Method,  proves  this  intuition  right.  For  all  but  a 
few  test  problems,  Anderson’s  method  is  clearly  superior  in  the  number  of  iterations 
required  for  convergence.  Even  for  the  exceptions,  Anderson’s  convergence  rate  was 
almost  identical  to  the  method  that  obtained  the  best  results.  In  addition,  there  was 
not  one  test  problem  for  which  Anderson’s  method  failed  to  converge  to  the  solution. 
There  were  test  cases,  Examples  3  and  8,  where  the  other  methods  had  the  problem 
of  limited  accuracy,  the  inability  to  achieve  convergence  with  a  precision  of  (715  even 
though  convergence  to  a  poorer  precision  is  obtained.  For  these  problems,  making 
variations  to  the  method  by  combining  different  relaxation  parameters  with  different 
amounts  of  discarded  iterates  still  did  not  achieve  6^15  convergence.  Therefore, 
it  is  the  author’s  conclusion  that  Anderson’s  method  will  consistently  solve  most 
numerical  problems  in  fewer  iterations  than  the  other  methods  studied  in  this  thesis, 
and  that  it  is  less  susceptible  to  limited  accuracy  than  are  the  other  methods. 

As  stated  previously,  it  should  be  emphasized  that  because  Anderson’s  method 
does  require  an  extrapolation  every  iteration  after  the  first  extrapolation,  the  com¬ 
puter  time  required  to  solve  some  fast  iterative  problems  may  be  more  than  if  the 
method  is  not  applied.  However,  for  most  problems,  especially  integral  equation 
problems  and  problems  with  a  divergent  generated  sequence,  fewer  iterations  is 
definitely  desired;  hence,  Anderson’s  method  will  usually  provide  the  best  results. 

Another  area  I  want  to  stress  is  the  reversing  of  the  generated  sequence  when 
applying  the  RRE  and  MPE  methods  to  a  convergent  generated  sequence.  Test 
results  show  that  this  procedure  will  produce  better  results  (though  there  were  a 
few  exceptions  for  certain  k  values  and  a  particular  problem)  than  if  the  sequence 
is  not  reversed.  For  a  divergent  sequence,  results  prove  that  the  original  sequence 
produces  the  best  result.  In  almost  all  test  problems,  the  RRE  and  MPE  methods 
gave  similar  results.  Even  though  these  two  methods  seldom  equaled  Anderson’s 
method,  they  consistently  outperformed  the  vector  Aitken  and  the  vector  e  methods. 
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There  axe  still  areas  of  study  that  can  be  investigated.  First,  the  GMRES 
method  of  Chapter  XII  can  be  fully  tested  and  compared  with  the  other  methods. 
A  second  area  that  can  receive  future  study  is  trying  to  convert  either  the  RRE  or 
MPE  methods  into  a  fully  dynamic  model.  By  using  the  principle  introduced  by 
Irons  and  Shrive  (1987)  in  Chapter  IV  for  the  scalar  case,  perhaps  a  fully  dynamic 
model  can  be  derived  for  the  RRE  and/or  the  MPE  methods. 
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APPENDIX 


CORRECTIONS  TO  ARTICLE  WRITTEN  BY 
SMITH,  FORD,  AND  SIDI  (1987) 

David  A.  Smith,  William  F.  Ford,  and  Avram  Sidi  wrote  an  article,  “Extrap¬ 
olation  Methods  for  Vector  Sequences,”  in  the  SIAM  Review,  Vol  29  (1987)  com¬ 
paring  acceleration  techniques.  These  methods  included  the  vector  epsilon  method 
with  both  types  of  inverses,  the  generalized  and  the  primitive;  the  MPE  method; 
and  the  RRE  method.  Several  of  the  comments  in  the  paper  concerning  their  test 
results  are  not  correct.  This  Appendix  details  the  errors  and  corrections  needed,  if 
appropriate. 

In  their  Example  2  (Example  2  in  Chapter  XI  also),  they  claim  that  the 
RRE  and  Vector  Epsilon  methods  “converge”  to  a  vector  approximately  equal  to 
(13.36,-1.940,5.532,-5.342).  This  is  not  the  case.  Both  methods  converge  to 
the  unique  solution  (1,1, 1,1).  When  converting  the  problem  to  the  Gauss-Seidel 
iteration  scheme  (88),  they  continued  to  use  the  original  vector  b  =  (10,4, 8, 6)r 
instead  of  the  converted  vector  ( D  +  L)~lb  =  (5, 1/3,  — 11/9, 163/9)r.  As  a  result 
they  determined  the  solution  of  a  different  fix  point  problem,  and  their  method 
converged  to  the  correct  solution  for  their  incorrect  problem.  Using  the  correct 
converted  vector,  the  RRE  and  Vector  Epsilon  methods  converge  nicely  to  the 
solution  (1, 1, 1, 1). 

They  also  state  that  because  the  system  has  a  zero  eigenvalue,  the  system  of 
equations  is  singular  for  k  =  4.  However,  for  the  initial  vector  they  used,  (0, 0, 0, 0), 
the  error  vector  does  not  lie  in  a  subspace  spanned  by  any  three  eigenvectors  of  the 
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iteration  matrix,  and  the  system  is  not  singular  for  this  starting  point  and  value 
of  k.  As  a  result,  the  MPE  method  is  exact,  in  the  absence  of  rounding  error,  for 
k  =  4  but  not  for  k  =  3.  If  the  initial  iterate  is  discarded,  as  discussed  in  Chapter 
XI,  then  k  =  3  is  appropriate. 

For  their  Examples  1  and  8  (Examples  3  and  4  in  Chapter  XI),  they  stated 
that  the  RRE  method  failed  to  converge  to  the  correct  solution.  Test  results  show 
that  the  RRE  method  does  converge  to  the  solution  in  both  cases.  For  the  first 
example,  convergence  is  very  similar  to  that  obtained  by  the  MPE  method  and  is 
obtained  for  all  values  of  k,  though  the  convergence  is  hampered  by  the  problem  of 
limited  accuracy.  The  second  example  is  a  problem  with  two  solutions.  The  RRE 
method  ( k  =  4)  converges  to  the  same  solution  as  the  basic  iteration,  (3,3, 3,3). 
However,  it  should  be  noted  that  for  k  <  4  the  RRE  method  caused  system  overflow 
for  this  problem. 
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